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Special Issue Information 


Dear Colleagues, 


Aircraft design is, as we know, the first fascinating step in the life of an aircraft, where visions 


are converted into reality. 


In a practical sense, aircraft design supplies the geometrical description of the aircraft. 
Traditionally, the output is a three-view drawing and a list of aircraft parameters. Today, the 
output may also be an electronic 3D model. In the case of civil aircraft, a fuselage cross- 


section and a cabin layout are provided in addition. 


In an abstract sense, aircraft design determines the design parameters to ensure that the 
requirements and constraints are met and design objectives are optimized. The fundamental 
requirements for civil aviation are payload and range. Many constraints come from 
certification rules demanding safety. The objectives are often of a financial nature, like 
lowest operating costs. Aircraft design always strives for the best compromise among 
conflicting issues. 


The design synthesis of an aircraft goes from the conceptual design to the detailed design. 
Frequently, expert knowledge is needed more than computing power. Typical work involves 
Statistics, the application of inverse methods, and use of optimization algorithms. Proposed 
designs are analyzed with respect to aerodynamics (drag), structure (mass), performance, 
stability and control, and aeroelasticity, to name just a few. A modern aircraft is a complex, 
computer-controlled combination of its structure, engines, and systems. Passengers demand 
high comfort at low fares, society demands environmentally friendly aircraft, and investors 


demand a profitable asset. 


Overall aircraft design (OAD) comprises all aircraft types in civil and military use, considers all 
major aircraft components (wing, fuselage, tail, undercarriage) as well as the integration of 
engines and systems. The aircraft is seen as part of the air transport system and beyond 
contributing to multimodal transport. Aircraft design applies the different aerospace 
sciences and considers the aircraft during its whole life cycle. Authors from all economic 
sectors (private, public, civic, and general public) can submit to this Special Issue (Sl). 


Education and training in aircraft design is considered as important as research in the field. 


The SI can be a home for those active in the European Workshop on Aircraft Design 
Education (EWADE) or the Symposium on Collaboration in Aircraft Design (SCAD), both 
independent activities under the CEAS Technical Committee Aircraft Design (TCAD). Please 
see http://AircraftDesign.org for details. 


Following the successful initial Special Issue on “Aircraft Design (SI-1/2017)”, this is already 
the second SI “Aircraft Design (SI-2/2020)”. Depending on the need, further special issues 
may follow. Activities in the past showed that aircraft design may be a field too small to 
justify its own (subscription-based) journal. A continuous open access special issue may fill 
the gap. As such, the Special Issue “Aircraft Design” can be a home for all those working in 


the field who regret the absence of an aircraft design journal. 


The Special Issue “Aircraft Design” is open to the full range of article types. It is a place to 
discuss the "hot topics" (zero-emission airplanes, electric flight, urban air mobility—you 
name it). The classic topics in aircraft design remain: 

e Innovative aircraft concepts 

e Methodologies and tools for aircraft design and optimization 

e Reference aircraft designs and case studies with data sets 
It is up to us as authors to shape the Special Issue “Aircraft Design” according to our interests 


through the manuscripts we submit. 


Prof. Dr. Dieter Scholz and Prof. em. Egbert Torenbeek, Guest Editors 
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Abstract: The article looks at publishing options in the field of aircraft design to find that no dedicated 
journal on aircraft design exists. For this reason, a Continuous Special Issue Aircraft Design of the 
well established journal “Aerospace” at the Open Access publisher MDPI is started. Often special 
issues of a journal are introduced for “hot topics”. Here, the subset “special issue” is used for a 
scientific domain—in this case “aircraft design”. Recurring single special issues are numbered in 
sequence and are identified by the year of the deadline for manuscript submissions. This allows for 
the delivery of several single special issues over time in a row without the need to define a publishing 
schedule up front. Together these single issues form the Continuous Special Issue Aircraft Design 
and offer a new publishing home for the aircraft design community. 


Keywords: aerospace; aviation; aeronautics; airplanes; aircraft; design; publishing; open access; 
special issue; MDPI; permalink; archive 


1. Introduction 


Aerospace consists of aeronautics (atmospheric flight) and astronautics (space flight) [1]. 
The scientific foundation of aeronautics is called aeronautical science, which includes aeronautical 
engineering. One of the many disciplines in aeronautical engineering is aircraft design. Aircraft design is 
the very first step in aeronautical engineering, where requirements are converted into a geometrical 
description of the aircraft. 

A list of scientific journals in the area of “aerospace” is available from Google Scholar [2], CWTS 
Journal Indicators [3], and the University of Illinois [4]. A ranking of the top 40 aerospace journals was 
done by the German Society for Aeronautics and Astronautics (DGLR) in 2015 [5]. Among all of these 
listed journals, only two carry “aircraft” in their name. These two are the Journal of Aircraft by AIAA 
(ISSN 0021-8669) and Aircraft Engineering and Aerospace Technology by Emerald (ISSN 1748-8842). None 
of them specifically deals with aircraft design. 

In 1998 the situation was much the same as today. For this reason, Prof. Egbert Torenbeek [6] 
started the journal “Aircraft Design” (ISSN 1369-8869) together with Prof. Dr. Jan Roskam [7] at 
Elsevier [8]. Both of them served as Editor. Torenbeek took care of authors in Europe and Roskam 
likewise took care of authors in the USA [9]. The journal started successfully and published 58 
articles in the four years until 2001 [10]. The subscription-based publishing model, however, proved 
inadequate to serve the rather small aircraft design community. As a consequence, Elsevier’s title had 
to be discontinued in 2002, when it was decided that the number of subscriptions was too low. 

Today, this pitfall can be avoided with the Open Access (OA) publishing model [11,12] because 
the papers finance themselves with article processing charges (APC), which is referred to as “gold OA”. 
In turn, the papers can be read free of charge on the Internet (“gratis OA”) and are in addition free 
of most copyright and licensing restrictions (“libre OA”). Often, the Creative Commons Attribution 
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(CC BY) license is used, for example, this article makes use of CC BY. Please refer to the bottom of this 
article for the CC BY logo and the link leading to further details. 

Furthermore, the motivation back in 1998 for establishing a journal was much the same as today. 
The European Workshop on Aircraft Design Education (EWADE) had been started in 1994 and was 
held every two years [13]. No workshop passed without expressing regret about the absence of a 
dedicated aircraft design journal [13-15]. In the same way as in the past, today people working in the 
domain of aircraft design also see each other as a close-knit and affectionate community. The European 
Workshop on Aircraft Design Education (EWADE) [16] was continued along with the Symposium on 
Collaboration in Aircraft Design (SCAD) [17], both became independent activities under the CEAS 
Technical Committee Aircraft Design (TCAD) [18]. CEAS is the Council of European Aerospace 
Societies [19]. Presently, CEAS has twelve European national aerospace societies as members. In 
eastern or central Europe engineers come together in a workshop called Research and Education in 
Aircraft Design (READ) [20]. Similarly, in the USA the AIAA Aircraft Design Technical Committee is 
very active [21]. AIAA is the American Institute of Aeronautics and Astronautics. For these and other 
aircraft design related communities a Special Issue Aircraft Design can serve as a publishing home in 
the absence of a dedicated journal for the discipline. 


2. Set Up of the Special Issue “Aircraft Design” 


MDPI (Multidisciplinary Digital Publishing Institute) uses the Open Access publishing model 
for all of its journals including its special issues. “Aerospace” (ISSN 2226-4310) is one such journal at 
MDPI which was started in 2014, and its aims include the “design... of aircraft” [22]. The first Special 
Issue Aircraft Design (SI-1/2017) [23] at “Aerospace” was managed by Guest Editor Dr. Mohammad 
Sadraey and includes five papers. The second Special Issue Aircraft Design (SI-2/2020) [24] is managed 
by Prof. Dr. Dieter Scholz, MSME [25] as Guest Editor. In addition, Prof. em. Egbert Torenbeek 
contributes with his experience as Honorary Guest Editor. 

As can be seen above, the special issues are consecutively numbered. The year given after the 
issue number is the year in which the Special Issue closed or is scheduled to close. This allows for 
the delivery of a continuous sequence of special issues over time without the need to define a certain 
publishing sequence from the start. It enables the Continuous Special Issue Aircraft Design, which 
stands for the sum of all individual Special Issues Aircraft Design and the mechanism to make this an 
ongoing activity. 

Each special issue has its own home page at MDPI. These separate home pages allow for changes 
to be expressed in the aims and scope or among the Guest Editors. In addition, the special issue home 
page can be seen as the front matter of a particular special issue. The banner of the special issue 
(Figures 1 and 2) can be shown on the home page as the graphical abstract of an Editorial located at the 
start of the list of papers. 
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Figure 1. Banner of the Special Issue Aircraft Design (SI-2/2020). 
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Figure 2. Banner of the Special Issue Aircraft Design (SI-1/2017). 


The URL of the currently open Special Issue Aircraft Design is always: https://www.mdpi.com/ 
journal/aerospace/special_issues/aircraft_design. Linking to this URL also means linking to the central 
anchor point of the Continuous Special Issue Aircraft Design as a whole now and in the future. This 
page shows links to direct readers to all past special issues. Linking directly to a single past issue is 
also possible: a long user-friendly URL is available as well as a short URL based on MDPI’s internal 
number of that special issue. Table 1 shows the systematic of these URLs. 


Table 1. Systematic of URLs of the Continuous Special Issues Aircraft Design. 


Closed SIs: Special Issue Aircraft Design (SI-1/2017) 

long URL https://www.mdpi.com/journal/aerospace/special_issues/aircraft_design_1_2017 
short URL https://www.mdpi.com/si/6497 2 

Open SI: Special Issue Aircraft Design (SI-2/2020) 

long URL https://www.mdpi.com/journal/aerospace/special_issues/aircraft_design ! 
short URL https://www.mdpi.com/si/25829 2 


' The currently open Special Issue Aircraft Design has always this URL. ? The internal special issue number can be 
obtained from the editorial office. 


Presently, the special issue is also advertised on a CEAS web page. The URL of the page is 
particularly easy to remember: http://journal.AircraftDesign.org. The page provides further details 
about the Continuous Special Issues Aircraft Design, e.g., related to possible reductions of the standard 
article processing charges (APC) [26] and links to all relevant web pages. 

The Digital Object Identifier (DOD) is always used when linking to an individual paper. The DOI for 
a recent “Aerospace” paper at MDPI is structured as follows: https://doi.org/10.3390/aerospaceVMMxxxx. 
Here, 10.3390 stands for MDPI, V is the Volume (related to the year, starting with 1 in 2014), MM is the 
month, and xxxx is the number of the paper. Allocation starts with 0001 in the beginning of each year, 
counting up. The setup is shown in Figure 3. 

When writing for the Special Issue Aircraft Design, authors should make sure to use persistent 
links to archived resources. Often we refer to journal articles that are usually long term archived by 
the publisher and have a persistent identifier (a DOI) to connect to the online resource. But in all other 
cases, it is the author’s responsibility to archive the web pages, PDFs, or other data files. Please refer to 
Appendix A to see how references can be archived with available tools on the Internet. 

A Continuous Special Issue as described above is not much different from what is otherwise 
known at MDPI as a “Topical Collection”. A topical collection is a variant of a special issue. Topical 
collections run continuously (without any numbering) under the same home page. Topical Collections 
are just another form to structure journal content. 
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SI-1/2017 


https://doi.org/10.3390/aerospace3040035 
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Figure 3. Identification of articles in the Continuous Special Issue Aircraft Design of the journal 
Aerospace by means of their DOI. 


3. Aims and Scope 


The Special Issue “Aircraft Design” is open to the full range of article types. In addition to original 
research articles, review papers, letters or communications, technical reports, and the extended version 
of conference papers are also accepted. Furthermore, an interest exists in aircraft design education. 
Certainly, the special issue is also the place to discuss topics like zero-emission airplanes, electric flight, 
urban air mobility—you name what is currently debated. Nevertheless, the classic topics in aircraft 
design remain: 


e Innovative aircraft concepts (Figure 4); 

e Methodologies and tools for aircraft design and optimization; 
e Reference aircraft designs and case studies with datasets; and 
e Aircraft design education. 
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Figure 4. Examples of innovative aircraft concepts [27]. 


The keywords are: aircraft, design, overall aircraft design (OAD), configuration, requirements, 
payload, range, certification, safety, constraints, objectives, synthesis, optimization, aerodynamics, 
drag, high-lift, structure, mass, performance, stability, control, aeroelasticity, engine, systems, operating 
costs, direct operating costs (DOC), passenger, cabin, ticket, price, environment, profit, asset, wing, 
fuselage, tail, undercarriage, landing gear, engine, systems. 

Authors from all economic sectors (private, public, civic, and general public) can submit to this 
Special Issue. Education and training in aircraft design is considered as important as research in 
the field. 


4. The Journal “Aerospace” at MDPI 


“Aerospace” (ISSN 2226-4310) is a well reputed journal as can be seen from the authors publishing 
with “Aerospace”. Its latest journal metrics CiteScore (CS), SRJ and SNIP (from Scopus, Elsevier) is 
given on the journal’s home page. Articles have a high visibility; papers are visible Open Access 
at the journal “Aerospace” and also alongside the Special Issue Aircraft Design as soon as they are 
ready. The journal “Aerospace” is covered by many databases [28] including Web of Science (Clarivate 
Analytics) and Scopus (Elsevier). Papers from “Aerospace” are archived for centuries to come at 
CLOCKSS and in the Swiss National Library Digital Archive [28]. “Aerospace” adheres to best practice 
in Open Access publishing (accessibility, openness, discoverability, reuse author rights, and many 
other criteria). This is expressed at the Directory of Open Access Journals (DOAJ) with the “DOAJ 
Seal” given to “Aerospace” [29,30]. MDPI is a member of many relevant publishing organizations 
(OASPA, COPE, STM, ... ) [31]. Membership to most organizations is only granted after a thorough 
check of the publisher and its journals. 

The journal “Aerospace” is known for rapid publication [32]. Manuscripts are peer-reviewed 
and a first decision is provided to authors approximately three weeks after submission; the length 
of the peer review itself can vary considerably, but reviewers are reminded by the editorial office to 
make the review a priority; acceptance to publication is undertaken in one week. Once accepted, 
the manuscripts undergo professional copy-editing, proofreading by the authors, final corrections, and 
publication on the journal website. This means that papers will be visible alongside with the “Special 
Issue Aircraft Design” and the journal “Aerospace” as soon as they are ready. 


5. Aircraft Design 


Aircraft design is the first fascinating step in the life of an aircraft, where visions are converted 
into reality. 

In a practical sense, aircraft design supplies the geometrical description of the aircraft. 
Traditionally, the output is a three-view drawing and a list of aircraft parameters. Today, the output 
may also be an electronic 3D model. In the case of civil aircraft, a fuselage cross-section and a cabin 
layout are provided in addition. 
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In an abstract sense, aircraft design determines the design parameters to ensure that the 
requirements and constraints are met and design objectives are optimized. The fundamental 
requirements for civil aviation are payload and range. Many constraints come from certification rules 
demanding safety. The objectives are often of a financial nature, like the lowest operating costs. Aircraft 
design always strives for the best compromise among conflicting issues. 

The design synthesis of an aircraft goes from conceptual design to detailed design. Frequently, 
expert knowledge is needed more than computing power. Typical work involves statistics, 
the application of inverse methods, and use of optimization algorithms. Proposed designs are 
analyzed with respect to aerodynamics (drag), structure (mass), performance, stability and control, and 
aeroelasticity, to name just a few. A modern aircraft is a complex, computer-controlled combination of 
its structure, engines, and systems. Passengers demand high comfort at low fares, society demands 
environmentally friendly aircraft, and investors demand a profitable asset. 

Overall aircraft design (OAD) comprises all aircraft types in civil and military use, considers all 
major aircraft components (wing, fuselage, tail, undercarriage) as well as the integration of engines and 
systems. The aircraft is seen as part of the air transport system and beyond contributing to multimodal 
transport. Aircraft design applies the different aerospace sciences and considers the aircraft during its 
whole life cycle [33]. 


6. Summary 


A journal titled “Aircraft Design” was published successfully from 1998 to 2001 by Elsevier, 
but had to be discontinued due to the low number of subscriptions. The demise was caused by a 
publishing model not adequate for the small aircraft design community. No other attempt to start 
an Aircraft Design journal has been made since then. The Open Access publishing model is a viable 
and better alternative for small communities. In a new endeavor toward creating something like an 
Aircraft Design journal, the subset of an Open Access journal within the wider topic “Aerospace” 
was used. The form of a Continuous Special Issue was chosen as the journal subset. The established 
journal “Aerospace” helps to overcome the problem of achieving “critical mass” for the new venture. 
Furthermore, MDPI provides a proven publishing infrastructure and support. 


Funding: This article received no external funding. 
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Abbreviations 

AIAA American Institute of Aeronautics and Astronautics (www.aiaa.org) 

CC Creative Commons (www.creativecommons.org) 

CEAS Council of European Aerospace Societies (www.ceas.org) 

CLOCKSS Controlled LOCKSS (www.clockss.org) 

COPE Committee on Publication Ethics (www.publicationethics.org) 

CS CiteScore [34] 

CWTS Centre for Science and Technology Studies (www.cwts.nl) 

DCLR Deutsche Gesellschaft fiir Luft- und Raumfahrt Lilienthal-Oberth e.V. 
(www.deglr.de) 

DOC Direct Operating Costs 

DOI Digital Object Identifier (www.doi.org) 

EWADE European VIOESHOE on Gonctare Design Education 
(www.ewade.aircraftdesign.org) 

ISSN International Standard Serial Number (www.issn.org) 

LOCKSS Lots of Copies Keep Stuff Safe (www.lockss.org) 

MDPI Multidisciplinary Digital Publishing Institute (www.mdpi.com) 

OAD Overall Aircraft Design 

OASPA Open Access Scholarly Publishers Association (www.oaspa.org) 

READ Research and Education in Aircraft Design (www.read.aircraftdesign.org) 
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SCAD Symposium on Collaboration in Aircraft Design (www.scad.aircraftdesign.org) 
SI Special Issue 

SNIP Source Normalized Impact per Paper [34] 

SRJ SCImago Journal Rank [34] 

STM Here: The worldwide association of STM publishers (www.stm-assoc.org); 
STM Science, Technology and Medicine 

TCAD CEAS Technical Committee Aircraft Design (www.aircraftdesign.org) 

URL Uniform Resource Locator 

Appendix A 


The Appendix explains the use of persistent links to archived resources as applied in this Editorial. Authors 
writing for the Special Issue Aircraft Design are encouraged to follow this practice! 

Journal articles (like those from MDPI) are usually archived and have a persistent identifier in the form of a 
DOI to connect to the online resource. But often, we make use of web pages, PDFs, or other data files taken from 
places on the Internet, where the resource is not archived and provided with a persistent identifier. We know from 
experience, “websites change, go away, and get taken down. When linked citations lead to broken, blank, altered, 
or even malicious pages, that’s called link rot.” [35]. Today, with tools at hand, it is the author’s responsibility 
to create a permanent, reliable, unbreakable link to an unalterable, archived record of the web page or the web 
resource cited in the work if this is otherwise not available. A comfortable way to do this is with Perma.cc 
(www.perma.cc). Perma.cc requires an account, which can be obtained free of charge e.g., from a participating 
university library. Alternatively, the Internet Archive with its Wayback Machine (www.web.archive.org) and the 
function “Save Page Now” can capture a web page or another online resource as it appears now for use as a 
trusted citation in the future. Links to the Internet Archive tend to be very long and would need to be shortened. 
This can be done e.g., with Bitly (www.bit.ly). The resulting short links fit well into the List of References. 

Following best practice, the List of References should have two links for each entry to a self-archived web 
resource; the original link and the link to its archived version. At MDPI typesetting rules only allow specifying 
one link. It is helpful that Perma.cc as well as the Internet Archive show also the link to the original resource in 
their archived documents. Therefore, specifying only the permanent link is not the best solution, but sufficient. 
If required, the reader would need to go via the archived version to find the original link. Please consult the 
References below to see how it works. 
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Abstract: Amphibious aircraft designers face challenges to improve takeoffs and landings on both 
water and land, with water-takeoffs being relatively more complex for analyses. Reducing the water- 
takeoff distance via the use of hydrofoils was a subject of interest in the 1970s, but the computational 
power to assess their designs was limited. A preliminary computational design framework is devel- 
oped to assess the performance and effectiveness of hydrofoils for amphibious aircraft applications, 
focusing on the water-takeoff performance. The design framework includes configuration selections 
and sizing methods for hydrofoils to fit within constraints from a flying-boat amphibious aircraft 
conceptual design for general aviation. The position, span, and incidence angle of the hydrofoil are 
optimized for minimum water-takeoff distance with consideration for the longitudinal stability of 
the aircraft. The analyses and optimizations are performed using water-takeoff simulations, which 
incorporate lift and drag forces with cavitation effects on the hydrofoil. Surrogate models are derived 
based on 2D computational fluid dynamics simulation results to approximate the force coefficients 
within the design space. The design procedure is evaluated in a case study involving a 10-seater 
amphibious aircraft, with results indicating that the addition of the hydrofoil achieves the purpose of 
reducing water-takeoff distance by reducing the hull resistance. 


Keywords: amphibious aircraft; hydrofoils; takeoff performance; computational fluid dynamics; 


optimization 


1. Introduction 


Amphibious aircraft have the potential to play an important role in passenger trans- 
port as part of general aviation, particularly in short-range flights [1]. Aircraft designed 
for takeoff and landings solely on runways are limited by the number of airports present 
in regions. The use of water bodies and ports as additional takeoff and landing points 
pose larger versatility and scope for missions with the use of amphibious aircraft. Mission 
profiles, weight estimation, fuel efficiency, range, payload, and stability are key considera- 
tions in the preliminary design and development of any aircraft. However, the takeoff and 
landing flight segments for amphibious aircraft present greater complications than ones 
in general aviation aircraft design because of the increased complexities involving water 
analyses. These aircraft suffer from several restrictions in analyses, due to the complexities 
of characteristic aerodynamic and hydrodynamic parameters, that make it difficult to 
non-dimensionalize and test models with respect to a length scale; these complexities are 
described in Section 2.1. 

Amphibious aircraft are mostly used as utility vehicles for search and rescue operations, 
payload delivery in remote and undeveloped areas, and firefighting efforts, among others. 
Sikorsky Aircraft developed numerous amphibious aircraft from the mid-1920s till 1940 for 
varied customers (Popular Aviation, October 1931: https:/ /books.google.com/books?id= 


https://www.mdpi.com/journal/ aerospace 


Aerospace 2021, 8, 10 2 of 31 


gZhdxqqWu4egC&pe=PA89). The Piaggio P. 136 was an amphibious aircraft used by the Ital- 
ian Air Force in the 1950s (http: / /www.aeroflight.co.uk/waf/italy /af/ital-af2-all-time.htm). 
More recently, the ShinMaywa US-2 has been in production, with a history of amphibious 
aircraft development since the 1950s (https:/ /www.shinmaywa.co.jp/aircraft/english/us2 
/us2_history.html). The Dornier Seastar is an amphibious aircraft from the 1970s with a 
new generation currently under development (https: //www.flightglobal.com/business- 
aviation /dorniers-new-generation-seastar-makes-maiden-sortie/137601.article). The AVIC 
AG-600 “Kunlong”, currently under development in China, is claimed to be the largest 
amphibious aircraft developed for aerial firefighting, maritime patrol, and search and res- 
cue operations (http: //www.xinhuanet.com/english/2018-05/13/c_137176292.htm). The 
research on amphibious aircraft development is lacking partly due to its focus on specific 
missions for utility purposes [2], as well as the increased complexities of hull and float 
analyses. The early developments of amphibious aircraft relied mainly on experimental 
water tank tests and empirical methods. As such, a thorough design exploration and opti- 
mization was not possible due to computational limitations [3,4]. A common strategy is to 
take conventional aircraft and convert them into amphibious configurations by the addition 
of a hull or float, which results in suboptimal designs. The aim for truly optimal amphibious 
aircraft design would require a systematic design procedure that considers preliminary 
sizing, stability and control, and shape optimization. A conceptual design and sizing frame- 
work for amphibious aircraft was developed by Cary [5]. While numerical optimization 
has been commonly used in aircraft design since the late 1970s [6,7], its application in am- 
phibious aircraft is still limited. Qiu and Song recently performed a response-surface-based 
multiobjective optimization to find the optimum hull step configuration, achieving an 18% 
improvement in takeoff distance without sacrificing cruise performance substantially [8]. 
Puorger et al. performed an aerostructural optimization for a fire-extinguisher amphibious 
aircraft [9], focusing mainly on the ground effect during low-altitude cruise, without any 
water-takeoff considerations. A water-takeoff performance calculation method using low- 
cost empirical models has also been recently developed by Wang et al. [10]; this method 
was derived based on digital virtual flight, and included a pilot model to study the impact 
of a lack of visual references during water-takeoff. The model was validated against some 
experimental data with less than 10% error. 

A particular technology for the potential improvement of amphibious aircraft water- 
takeoff performance is the hydrofoil. In this paper, a hydrofoil is defined as a lifting surface 
that travels through water. The implementation of the technology allows a ship’s hull to 
reach a planing stage more quickly, and effectively unport the hull from water. This reduces 
motor effort by reducing hull resistance, which allows it to travel at higher speeds. They 
have been researched and implemented in water-based craft since the late 1800s, with exten- 
sive research projects performed between the 1930s and the 1960s to improve performance 
of marine vehicles [11]. Experiments showing the effectiveness of hydrofoils dated as early 
as 1861 by Thomas Moy [12], in 1907 and 1914 by the Wright brothers [3], and in 1937-1940 
by Sottorf [13], among many others. Most modern studies of hydrofoils relate to marine 
and naval technologies, including unmanned surface vehicles (USV) [14-16], aquatic un- 
manned aerial vehicle (AquaUAV) [17], turbomachinery [18], naval propellers [19], the 
propulsion of marine vehicles [20,21], or axial-flow pumps [22]. Hydrofoil boats, such as 
the “Raketa” developed in Russia in 1957, were manufactured and used for commercial 
transporation [23]. Studies of hydro-aerodynamic characteristics of hydrofoil systems in 
vessels were performed by Plisov and Rozhdestvenski [24]. A hydrofoil, as opposed to 
an airfoil, can suffer from cavitation and ventilation [25], further described in Section 2.2. 
Shen and Eppler researched various hydrofoil profiles using inverse design methods for 
U.S. hydrofoil craft [26]. Studies of two-phase boundary layer control in cavitating and 
supercavitating conditions were performed by La Roche and Trevisani in relation to the 
Supramar company [27]. Garg et al. performed multipoint hydrodynamic shape and hy- 
drostructural optimizations of 3D hydrofoils in partially cavitating conditions to improve 
fuel efficiency of ships [28]. This analysis, however, did not consider multiphase flows 
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and enforced a cavitation constraint based on the pressure coefficient to obtain a shape 
that did not undergo cavitation, which might not be suitable for high-speed conditions. 
Vernengo et al. performed shape optimizations of 2D supercavitating hydrofoils with 
multiphase flows for the development of high-speed planing craft [29]. 

Amphibious aircraft with hulls can also benefit from hydrofoils similarly by reaching 
their takeoff speeds more quickly, reducing the distance and time in the water-takeoff 
stage. The research on hydrofoils for amphibious aircraft design is sparse; the number 
of designed hydrofoil profiles is far lesser compared to that of airfoils, and the number 
of investigations into hydrofoil shape optimization with the relevant constraints is much 
lesser, as well. The Thurston Aircraft Corporation published a report which extensively 
considered hydrofoils for seaplane design [3] based on experimental results, but not for 
amphibious aircraft with the relevant constraints. David Thurston also investigated general 
amphibious aircraft design with minor focus on the implementation of a hydrofoil in a 
working concept called the HRV-1 flown in 1974, but the rigor of analysis for the takeoff 
condition was limited to a basic sizing and did not consider multiple configurations [2]. 
Cary applied Thurston’s procedure into the amphibious aircraft design framework with 
hydrofoil considerations, but the scope of foil profile selection, viability of design, and 
stability were not investigated [5]. The hydrodynamic lift for flying boats and seaplanes 
via the use of hydrofoils for a particular retractable configuration was invented by Zimmer 
and patented under Dornier Luftfahrt [30]. The LISA Akoya (LISA Airplanes: http:// 
lisa-airplanes.com/en/light-amphibious-aircraft-akoya/), which is a two-seater aircraft 
designed for leisure flight, is an amphibious aircraft known to implement this technology 
in its design. 

A rigorous computational framework that can systematically evaluate and assess the 
impact of hydrofoil performance in amphibious aircraft is still lacking. The numerical 
studies on amphibious aircraft performance mentioned above [8-10] did not consider 
hydrofoils and thus would not be directly applicable to amphibious aircraft with hydro- 
foils. Such a framework will enable tradeoff studies and optimization, which can help 
bring revolutionary technological improvement to amphibious aircraft design, e.g., by 
optimizing the shape of the hull and hydrofoils. The current investigation aims to formu- 
late a computational design framework for preliminary hydrofoil design in the context 
of amphibious aircraft application. This will include procedures for hydrofoil sizing and 
placement, with longitudinal stability considerations. One goal of the implementation is to 
provide a “riding” surface for amphibious aircraft while in the planing stage to reduce hull 
drag by minimizing its water contact. In particular, we will investigate and compare the 
performance of amphibious aircraft hull with and without hydrofoils in different speed 
regimes during the water-takeoff procedure. The results will give us physical insights into 
the design problem. 

The paper is structured as follows. Section 2 discusses the framework of fluid dynam- 
ics important in amphibious aircraft design, including hydrofoils. Section 3 outlines the 
research methodology, including the sizing procedure of the hydrofoil, the explanation of 
the water-takeoff model for amphibious aircraft, the considerations of the effects of the 
hydrofoil on the stability of the aircraft, and the design framework for the analysis and 
optimization of the hydrofoil. Section 4 discusses the application of the hydrofoil design 
framework to a 10-seater amphibious aircraft design with its results. The summary of key 
findings and conclusion of this work are then presented in Section 5. 


2. Theory and Hydrofoil Design Challenges 
2.1. Hydrodynamic Analyses 


Three key non-dimensional parameters in water analysis are the Reynolds (Re), Froude 
(Fr), and Weber (We) numbers: 








Re = Oe, Fr= : We = © ; (1) 
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where p is the density of the fluid, u is the freestream speed, L is the reference length of 
the body of analysis, p is the dynamic viscosity, g is the gravitational acceleration, and 
« is the surface tension between the two fluids under consideration. Fluid parameters 
corresponding to air will be denoted by a subscripted A, e.g., wa, and similarly W for 
water. The Weber number does not scale accordingly with Reynolds and Froude numbers 
because it goes as the square of the velocity, so the benefits of non-dimensional analysis are 
lost when attempting to size a water-based component of a transport vehicle based on a 
scale model with experimental results. The case worsens with the Froude number going as 
the reciprocal of the square root of the length against the Reynolds numbers’ proportional 
relationship to length. The mitigating strategy is to use the same Froude number in the 
model tests and to adjust for different Reynolds numbers when scaling. Some errors exist in 
water spray, wave pattern and foaming predictions due to the difference in Weber numbers, 
but these are negligible in resistance prediction of scaled-up hulls [31]. 

A preliminary calculation under assumptions of constant viscosity and density shows 
that the Reynolds number in fresh water is approximately 16 times greater than the 
Reynolds number in air for the same speed and reference length. Lift and drag forces 
in fresh water are approximately 815 times (ow/p, ~ 840 at International Standard At- 
mosphere sea level conditions with salinity considerations) greater for foils of the same 
profile and dimensions traveling at a given speed with completely attached flow in both 
cases. Turbulence, cavitation, and ventilation effects correspondingly result in non-trivially 
determined differences between the two mediums. 

Hydrodynamic forces important in hull design, such as resistance and buoyant force, 
are non-dimensionalized by division with pw¢B°, where B is the hull width, as opposed to 
the product of the dynamic pressure and wing area in aerodynamics; this is usually because 
the hull lengths are fixed and the widths are varied for design analyses. Archimedes’ 
principle also justifies the inclusion of gravitational acceleration as a term to account for 
buoyancy forces. 

The following hydrodynamic coefficients, called the load, resistance, and velocity 
coefficients (Ca, Cr, and Cy, respectively), play an important role in the water-takeoff 
analysis of an amphibious aircraft with a hull as non-dimensional measures of buoyancy A, 
resistance R, and speed u: 


A R Uu 
Ch =———, Crpe=— KH Cy= 
PwgB° pweB° /oB 


The trim angle trim of a hull is defined as the angle at which a boat must be longitudi- 
nally inclined for a given speed with respect to its orientation at rest such that it maintains 
optimal performance via minimization of resistance generated. The resistance and trim 
angle of a hull are experimentally determined with variations against speed via tank tests 
of models, which are then scaled using methods, such as one from the International Towing 
Tank Conference in 1978 (ITTC-78) [31]. 

The non-dimensional lift C,, drag Cp, and moment Cy coefficients that characterize 
airfoil performance are also used for hydrofoils with the freestream density and dynamic 
viscosity of water as reference. However, phase changes can occur over lifting bodies 
traveling in water at speeds within the takeoff speed regime of aircraft, leading to additional 
complications not seen in airfoils. Here, we use the generic non-dimensionalization for the 
lift L, drag D, and moment M: 





(2) 


c D M 


Cr = ——, Cp= 7 Cu=7——, 3 

:  ou2S © 5pu2S " 5pUu2Sc ©) 
where 0, u, S, c refer to density, velocity, area, and chord of the lifting component, respec- 
tively. The subscripts hf, h, and w will denote the parameters corresponding to hydrofoil, 
horizontal tail, and wing, respectively. 
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2.2. Cavitation and Ventilation 


Cavitation and ventilation can notably reduce hydrodynamic efficiency of hydrofoils. 
Cavitation occurs when the local pressure approaches the saturated vapor pressure [28], 
which can cause flow-induced noise and vibration, decreased lift and thrust, and increased 
drag [32]. Ventilation refers to the phenomenon where air bubbles are found in the sub- 
merged part of the body [15], and this can result in a loss of lift [3,16]. The measure 
of ventilation effects is considered to be the depth of submergence of the lifting surface 
from the free surface. In the case of a ventilation analysis for a hydrofoil, the depth of 
submergence in water from sea level is the appropriate measure. Early work to optimize 
the hydrodynamics of hydrofoils did not consider cavitation [33], though most subsequent 
works considered it [18,20,34-36]. These works mainly used potential flow as their bases. 
In recent decades, some studies have used high-fidelity models, in particular the Unsteady 
Reynolds-averaged Navier-Stokes (URANS) approach with a physical model of cavita- 
tion [19,37]. While the cavitation effects of hydrofoils have been well-investigated, studies 
of ventilation are still scarce [14,38]. As ventilation is a highly complex phenomenon for 
which systematic data are very difficult to obtain either experimentally or numerically, it is 
not explicitly considered in the preliminary design phase discussed in this paper. 

Cavitation is an unsteady phenomenon that takes place in water when the local 
pressure on a surface is below the saturated vapor pressure of water, with bubble and 
vapor formation taking place along the lifting surface. This is known to cause “cavitation 
damage” in the form of corrosion of rotor blades of boat motors. In the case of aircraft, a 
disadvantageous effect is the increase in drag of the aircraft during takeoff and possibly 
some form of cavitation damage along the hull for certain configurations. Cavitation 
performance is usually compared via the non-dimensional cavitation number, defined as: 


ey en ae A 


1 7 Poo = Patm + owgh, (4) 
TPWUoo 


where fou is the freestream pressure, Uoo is the freestream speed, Patm is the atmospheric 
pressure, py is the vapor pressure, and h is the reference height corresponding to the 
hydrostatic pressure head. 

The condition for cavitation inception over hydrofoils is given as —Cp,,,, = Ca, where 
Comin 1S the minimum pressure coefficient over the hydrofoil. As a result, cavitation 
inception usually occurs on the upper surface of foils as the local pressure is lower than the 
freestream pressure for foils under lifting conditions. Numerical studies in hydrodynamic 
analysis and optimization for foil shapes usually consider variations of lift and drag forces 
against cavitation numbers. 

Another issue involving cavitation inception is that the Reynolds, Mach, and cavitation 
numbers now play roles in the similarity analysis of flows over hydrofoils due to the flow’s 
additional dependence on the vapor pressure of the fluid in liquid state. If a particular 
Reynolds number is set for analysis, its cavitation and Mach numbers are determined for a 
fixed reference length. As the cavitation number is independent of a length scale, adjusting 
the reference length to obtain a freestream speed for a fixed Reynolds number adjusts the 
cavitation number accordingly. This changes the characteristics of cavitation formation 
over the hydrofoil, as the pressure coefficient conditions for cavitation inception change, 
viz. the location of incipient cavitation over the hydrofoil surface will change. 

The major issue with existing studies of hydrofoils is that they have originally been 
designed for cruising speeds of ships, which correspond to certain cavitation numbers as 
their design points. In the water-takeoff regime for an amphibious aircraft, such a fixed 
design point does not apply as the flow is accelerating. 


3. Research Methodology 


The work presented in this paper focuses on the hydrofoil design and analysis. The 
general preliminary design, weight estimation, and sizing procedures for the aircraft 
are therefore outside the scope of this work. For the discussion presented below, we 
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assume that the aforementioned procedures have been completed and the maximum 
takeoff weight, the required wing area, aircraft stall speed, and the powerplant selection 
have been determined. A hull selection and design are also already determined based on 
hull resistance and trim angle variations with speed. The Douglas Sea Scale [39] is used 
to measure the height and swell of a sea on a 0-10 degree scale. In this paper, a level 0 
(no wave) condition in calm water is assumed in all considerations for the initial design. 
The high-fidelity computational analyses are performed in 2D with a fixed reference length 
based on the sizing procedure. The 2D analyses are deemed suitable for the preliminary 
design stage [29]. Performing high-fidelity 3D analyses that involve multiphase flows could 
be computationally prohibitive and would be more appropriate at detailed design stages. 

The drag force caused by the empennage is considered as negligible in all equations 
of motion compared to the forces generated by the wing and hydrofoil. All aerodynamic 
forces are assumed to be concentrated as point loads from the aerodynamic center of the 
respective lifting surface, assumed to be located at 25% of the chord length from the leading 
edge. No control surfaces for the hydrofoils will be considered in the current investigation 
as the design framework considers the stability without the need of these devices. 

Based on these assumptions for the design scope, the hydrofoil conceptual design 
and preliminary sizing, water-takeoff analysis, and stability analysis are discussed in 
this section. 


3.1. Hydrofoil Conceptual Design Framework and Preliminary Sizing 


The following design procedure considers hull-based amphibious aircraft instead 
of float-based, as the primary goal is to minimize large hull contact with water. This 
section provides guidelines to determine the suitable hydrofoil configuration, geometry, 
and location for preliminary design purposes. 


3.1.1. Configuration and Profile Selection 


There are many possible hydrofoil configurations for amphibious aircraft including, 
but not limited to, surface-piercing hydrofoils, foldable-protruding hydrofoils, and strut- 
based hydrofoils [3,24]. As the name suggests, a surface-piercing hydrofoil has a hydrofoil 
protruding directly from the hull of the aircraft that pierces the water surface, such as the 
ones in LISA Akoya. When such a hydrofoil is foldable, it is considered as a foldable- 
protruding hydrofoil. The last category refers to a hydrofoil with a strut-based structure 
extending from the body of the aircraft. The surface-piercing and strut-based hydrofoil 
configurations are depicted in Figure 1. For the water-takeoff consideration in our design 
and analysis discussion, the hydrofoils are assumed to be fully deployed. The major 
dimensional constraints that will be discussed are relevant to this state, as well. 


Top view fo | 
PO 


Top view 


Waterline 


Waterline 







Hydrofoil 










Strut 
oO Hydrofoil 





ye 


Front view 
(a) (b) 
Figure 1. Hydrofoil configurations. (a) Surface-piercing Hydrofoil; (b) Strut-based Hydrofoil. 


Front view 


The number of hydrofoils on the aircraft is an important consideration in the weight 
and stability estimation of the initial design. The LISA Akoya implements two hydro- 
foils, one set near its wing and one set below the empennage; this configuration seems to 
eliminate the complex hull design considerations, such as stepped and planing configura- 
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tions, and enables the design of more aerodynamic fuselages as the hydrofoils provide lift 
throughout the water-takeoff run and reduce the wetted hull area. 

To select an appropriate hydrofoil profile, we identify several criteria to ensure efficient 
performance. These criteria include: low hydrofoil drag coefficient within the entire takeoff 
speed range, and high lift-to-drag ratios within the speed range in which hull resistance is 
dominant. To ensure stability, the rate of change of the hydrofoil moment coefficient with 
respect to angle of attack should be negative within the range of hull trim angle. A higher 
hydrofoil lift coefficient is desirable as it will reduce the required hydrofoil area, to be 
further discussed in Section 3.1.2. 

There are three main classifications for cavitation of hydrofoils, namely (1) subcavi- 
tating, where flow is fully attached over the hydrofoil; (2) partially cavitating, where flow 
separation of water takes place at some transition point with vapor cavities partially form- 
ing over the upper surface of the hydrofoil; and (3) supercavitating, where the hydrofoil 
undergoes separated flow from water with the formation of a large vapor cavity over 
the entire upper surface [3]. Past literature has shown hydrofoils have been designed to 
be efficient in either subcavitating or supercavitating regimes. Note that the assumption 
of the “aerodynamic center” of the hydrofoil located at approximately 25% of the chord 
length does not apply for these foils under cavitating conditions. The available databases 
of hydrofoils (e.g., University of Illinois at Urbana-Champaign (UIUC) Airfoil Database: 
https: / /m-selig.ae.illinois.edu/ads/coord_database.html) pale in comparison to their aero- 
dynamic counterparts and a designer must resort to sifting through historical literature 
to find shapes. Moreover, the cavitation phenomenon will change the flow performance 
characteristics in hydrofoils compared to those of airfoils, as will be seen in the lift and 
drag coefficient profiles presented in Section 4. 

The estimated range for the minimum pressure coefficient to ensure a subcavitat- 
ing regime for an amphibious aircraft can be determined by its design water-takeoff 
speed. Consider, for instance, the LISA Akoya, a small two-seater aircraft, which has a 
stall speed of us = 23 m/s and an estimated average hydrofoil height of 1m from the 
free surface (http: //lisa-airplanes.com/en/light-amphibious-aircraft-akoya/equipment- 
performance/). Estimating its takeoff speed at sea level by uto = 1.2, us = 27.6 m/s, its 
cavitation number constraint is the following: 


101,325 + 998.2 x 9.81 x 1 — 3640 
Be at si cia SAE OY) (5) 
VX 998.2 x 27.62 


This means that the lowest pressure coefficient over the hydrofoil must be Cy... = 
—0.2827 to prevent incipient cavitation; it is unlikely for subcavitating foils to satisfy this 
pressure coefficient constraint with beneficial (L/D) ratios. Thus, at the required speeds 
for water-takeoff, which are bound to be higher for larger aircraft, cavitation is almost 
an inevitable phenomenon as the cavitation number reduces with increase in dynamic 
pressure underwater, indicating that supercavitating hydrofoils are the only viable solution. 
Thurston and Vagianos showed that hydrofoils designed for supercavitating regimes 
far outperformed their subcavitating counterparts in terms of (L/D) ratios as a result 
of negligible effects of ventilation in supercavitating conditions, and should be primary 
considerations for amphibious aircraft [3]. A study on the performance of hydrofoils 
designed for subcavitating and partially cavitating regimes for use in amphibians supports 
this conclusion [1]. It is thus imperative to use supercavitating hydrofoils in amphibious 
aircraft applications. 

Experimental data exist for some supercavitating profiles [40], which can be used 
as a reference when sizing the hydrofoil. The attainable maximum lift coefficient for 
supercavitating hydrofoils appears to be Cy, © 0.27 including ventilation effects in the 
supercavitating regime according to Thurston’s review in the 1970s [2]. However, some 
supercavitating profiles, such as the SCSB profiles developed by Brizzolara [14], have 
shown improved performance in subcavitating and partially cavitating conditions in terms 
of (L/D)h z ratios compared to previous designs. 
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3.1.2. Preliminary Sizing 


The design considerations and criteria to determine hydrofoil geometry parameters in- 
cluding area, sweep angle, taper, dihedral, and aspect ratio are summarized and presented 
here. The selection is mainly based on the physical properties and hydrofoil performance. 
The span and incidence angle selection process, which relies on a water-takeoff minimiza- 
tion procedure, will be described separately in Section 3.1.3. 


3.1.2.1. Hydrofoil Area 


The hydrofoil needs to generate lift to unport the hull from the water while it is 
submerged. The stall speed of the hydrofoil (uj,,) is defined as the lift required by the 
hydrofoil to counter the weight of the aircraft minus the lift of the wing at the same speed. 
The required area for this condition can be derived from the results of the preliminary 
wing and horizontal tail sizing procedures by equating this stall speed to the speed at 
which the horizontal tail’s elevator (uy, ) is able to provide effective corrections as an initial 
design point. 


Unf, = Un,- (6) 


The required area 5j;,,, 1s then determined by equating the lift generated by the 
hydrofoil at its stall speed to the difference between the weight of the aircraft W and the 
lift generated by the wing in the same condition: 


1 1 
s Pw. Shtfreq Lit fimax — W - 5 PAUL, SwCLy (7) 


2W — Paus SwCL,, 
—— Sh freq — ee eee: ees 


(8) 


Z 
PW ie. Ch, Pas 


where W is the aircraft weight. A high C;,, fmax will therefore reduce the effective hydrofoil 
area, which reduces the drag generated by the hydrofoil. Note that this is the required area 
for a hydrofoil with no dihedral to generate sufficient lift for the aircraft weight, which can 
be slightly generalized as will be discussed in the dihedral sizing description below. 

The determination of C;,, depends on the Reynolds number, set by uy,, and ay, 
determined by the wing incidence angle ay, and hull trim angle at the corresponding point 
during the water-takeoff. As a first estimation in the design stage, C,,, can be determined 
via quick, low-fidelity analysis tools, such as XFOIL [41]. The value of C;,,,_ depends 
on the hydrofoil profile selection, as mentioned in Section 3.1.1. This method implicitly 
accounts for the losses of lift due to ventilation effects as a heuristic, assuming that the 
hydrofoil will not yield the ideal maximum lift coefficient at this design point due to 
ventilation. 


3.1.2.2. Sweep Angle and Taper 


King and Land [42] observed that the implementation of sweep in subcavitating 
conditions delayed the onset of cavitation inception to higher speeds for angles between 
0-—45° at the cost of decreased (L/D) ratios, so its implementation might be beneficial 
for aircraft designed for water-takeoffs in only subcavitating conditions for hydrofoils. It 
is, however, not efficient for supercavitating hydrofoils as they are designed for optimal 
performance in supercavitating conditions. They also observed that the implementation of 
taper has little effect on the characteristics of the flow and the performance. 


3.1.2.3. Dihedral Angle 


For surface-piercing hydrofoils, anhedral would be required to maintain lifting perfor- 
mance while the hull is unported from the water. The effective area of the hydrofoil with 
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a dihedral angle oy is then given by Sy¢ = Sy an / cos? 6), f, Which is then substituted into 
Equation (8) to give: 

2W — PAuy, SwCLy 
© PWH EE China, C8” Sif 


Shf (9) 

For a hull-protruding hydrofoil with a generic foil profile to be a beneficial lifting 
surface while reducing hull contact area, 0,¢ should be between —45° and 0° to enable 
unporting the hull while still being submerged, as the hydrofoil should protrude below the 
base of the hull and not above. 


3.1.2.4. Aspect Ratio 


The aspect ratio (AR) must be set, for a hull-protruding hydrofoil with positive 0), 
such that the height of the deployed hydrofoil is not greater than the distance of the landing 
gear ’s wheels from the aircraft’s center of gravity (CG, Z7¢). 


2 
A216 


ARns < O<a<l, (10) 


Sh freq tan? Onf ’ 


where a is a factor of safety that can be set by the designer to account for ground strikes 
during takeoffs on land in case of mechanical failures. 

The aspect ratio selected should ideally be high (i.e., greater than 5) according to 
Thurston [3]. For aircraft with strut-based hydrofoils, this may result in the effective span 
of the hydrofoil being larger than the fuselage width. In these cases, retractable/folding 
mechanisms should be considered in the design of the strut to fold the hydrofoil and retract 
it into the fuselage/hull. The determination of the aspect ratio is further investigated by 
optimization of the span length of the hydrofoil, discussed in Section 3.1.3. 


3.1.3. Span and Incidence Angle Optimization 


In this work, the hydrofoil’s span and incidence angle are determined via the min- 
imization of the water-takeoff distance, in which they are treated as design variables. 
The water-takeoff distance calculation, as the objective function, will be elaborated in 
Section 3.2. A lift constraint is considered, to model that the total lift provided by the 
lifting surfaces (wing, hydrofoil and horizontal tail) does not exceed the aircraft weight 
during the water-takeoff procedure, to satisfy the inequality constraint in residual form 
R(Ly¢) = Lag + Lw + Ly — W < 0. The variable bounds for the speed and angle of attack, 
which are used in the water-takeoff calculation, also need to be specified. The design 


limit of aspect ratio, which is related to the span as byr¢ = ,/ARyz x Spyz, is also imposed 


following the discussion in Section 3.1.2. The optimization procedure is summarized in 
Table 1. 


Table 1. Optimization problem formulation for the optimal hydrofoil span and incidence angle. 


Optimization 
Minimize 


Design variables 


Constraints 


Bounds 


Function Variables Description 
xw Water-takeoff distance. 
bn gf Span length of the hydrofoil. 
On f, Incidence angle of the hydrofoil. 
Lip + Lw + Ly <W The total lift provided by the lifting surfaces 
must not exceed the aircraft weight. 
O<u<urTO Speed range within takeoff regime. 
Xmin S &nf < &max Angle of attack within operational bounds. 


ARng_, < ARng < ARn Fann Aspect ratio within design bounds. 
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The corresponding extended design structure matrix (XDSM) (https://github.com/ 
mdolab/pyXDSM) representation [43] is shown in Figure 2. In this representation, math- 
ematical functions with well-defined inputs and outputs are colored in green, implicit 
components that solve residual equations are colored in red, optimizations and design 
of experiments are colored in blue, solvers pertaining to physics are colored in orange, 
and metamodels are colored in yellow. Data links between components are presented as 
thick, gray lines with slanted boxes denoting the communicated variables; inputs into 
components are fed “right-down” and outputs from components are returned “left-up”, 
depending on the order of execution. Inputs and outputs that are not connected to any 
other systems, such as initial guesses (denoted with —o subscripts) and optimization results 
(denoted with —* superscripts) of some optimizations, are shown as parameters inside 
white boxes. Processes are indicated with black lines representing arrows. 


W.Buotror fp [W] 


Takeoff Distance 







Minimizer 


Water- Takeoff 
Simulator 


Lift Constraint 


Figure 2. XDSM diagram for the hydrofoil span and incidence angle optimization. 


In this study, we only consider hydrofoils with rectangular planforms, as sweep and 
taper have no appreciable effects on the hydrofoil performance from the discussion in 
Section 3.1.2.2. To start the optimization procedure, an initial span length is selected based 
on the initial guess of the aspect ratio. This span length is then used to determine the chord 
length cy, of the hydrofoil that is used in computational fluid dynamics (CFD) analyses, 
elaborated in Section 3.4.2. 

The hydrofoil must be able to operate at the optimum angle for efficient performance 
in the takeoff regime. The angle of attack of the aircraft is subject to changes during 
water-takeoff due to hull designs and wave formations, which results in changes of the 
angle of attack of the hydrofoil. Hence, the incidence angle should be determined by 
accounting for the hull trim angle while considering the required angle of attack for the 
chosen hydrofoil profile: 


Kh f, = Crim — Xnf- (11) 


The angle of attack varies with time as a result of the hull’s dynamic trim angle. The 
consideration of an actuated system for adjusting the angle of attack, while possible, is too 
complex for an initial design consideration and is deemed beyond the scope of the present 
study. Note that such an implementation would add structural weight, which needs to 
be considered in the design process. The lift constraint is considered for the reason that if 
the aircraft and the hydrofoil are unported from water, the loss of lift from the hydrofoil 


11 of 31 


Aerospace 2021, 8, 10 


is extremely large due to the change in density of the fluid from water to air. This would 
result in the weight exceeding the total lift in this condition, and the aircraft would crash 
onto the water, resulting in unstable effects on the resistance and trim. 


3.1.4. Location Optimization 


The hydrofoil location in the amphibious aircraft is determined primarily based on the 
stability and trim requirements. The hydrofoil must generate minimal moments about the 
aircraft center of gravity (CG) until the speed at which the elevator is effective enough to 
provide corrections is reached, while simultaneously providing optimal lift to minimize hull 
resistance. This objective is attained by considering minimization of the maximum stabilizer 
force required during a simulated water-takeoff run. The reasoning for the consideration 
of this objective function will be presented in Section 3.3, and the water-takeoff simulation 
model will be described in Section 3.2. The optimization problem formulation to determine 
the hydrofoil location is summarized in Table 2. The hydrofoil location is indicated by the 
coordinates of its reference point, (xj,¢, Zn¢) with respect to the CG as origin. The optimum 
span and incidence angle from the takeoff distance optimization problem in Section 3.1.3 
are provided to the stabilizer force minimization problem, which determines the optimal 
hydrofoil position. The XDSM diagram is presented in Figure 3. 


Table 2. Optimization problem formulation to determine the hydrofoil’s location. 


Optimization 
Minimize 
Design variables 


Bounds 









(Uh ps Zp) 


nf, On 


LWw Lal, er /] 


Stabilizer Force 


Function Variables Description 
max|Ly| Maximum horizontal stabilizer force required during takeoff. 
(Xnp, Znf) Coordinates of hydrofoil’s reference point. 
Xnose < Xnf < 0 Horizontal position between aircraft nose and CG. 
Z1g S Zn < Zpase Vertical position between landing gear and aircraft base. 


Gn Timp OF 
(nf, 2nf) 





Minimizer 


Takeoff Distance 






Minimizer 
Water- Takeoff 
Sanne 
Simulator 


Figure 3. XDSM diagram for the hydrofoil location optimization. 


3.2. Water-Takeoff Analysis 


In this analysis, the water-takeoff distance is defined as the horizontal distance traveled 
by the aircraft during acceleration in water, to reach a pre-determined speed at which the 
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aircraft generates sufficient lift to achieve takeoff at a particular angle of attack. The 
water-takeoff analysis is performed based on Newton’s second law, as follows: 


mu = T cos ¢, (12a) 
ee 
Thrust 
_ PwsB Ch Cha Chinn 





Ca Ca 
— BS="°.Cp, where —> =1—‘ -°2— : 12b 
ie Oe Ca, 2(W/Sw) Cay Ce 


Wave resistance 


1 C 
— pPw Swell? Cr, (12c) 
~ 
Viscous resistance 


uz 


we} 


~ 


pa(CD, COS Ay + Cy, SiNaw) + Pw (Co,, COS Ay ¢ + Cy, Sin ang) | F (12d) 





Contributions from wing and hydrofoil lift and drag 


where m and u denote mass and acceleration, respectively. Equation (12d) consists of 
the contributions towards lift and drag from the wing and hydrofoil at their respective 
angles of attack, and Equation (12a) is the thrust force T at some thrust inclination angle ¢;. 
The remaining force terms on the right-hand side, such as the hull resistance and viscous 
resistance, are described in detail below. 

Equation (12b) represents the wave resistance (unrelated to wave drag) of the hull 
with respect to variations in its trim angle based on towing tank tests. This term is scaled 
by the submerged volume ratio Ca/C,, at a given instant during takeoff with respect to 
the submerged volume at rest. At each point of the water-takeoff calculation, the ‘effective 
weight’ of the aircraft (its weight minus the total lift) is considered to be balanced by the 
buoyant forces, hence A = W — > Lj, 1 € {hf,h,w}, which is non-dimensionalized to give 


1 

Ca. During the planing stage, the hull contact with water is minimal, i.e., Ca/Ca, © 0, 
where Ca, is the maximum load coefficient in the absence of lift, but it may generate 
viscous skin-friction resistance proportional to a small wetted area Swet in contact with the 
water, hence an additional viscous resistance term is added for the planing stage given 
by Equation (12c). Determining this wetted area is complex without knowledge of the 
exact waterline height and the spray pattern, so an initial value is set based on a geometric 
estimation and then scaled by the submerged volume ratio Ca /Ca,. The skin-frictional 
resistance coefficient Cr is obtained via the ITTC-57 formula, based on data from empirical 
towing tank tests [31]: 


07 
Ce= ee (13) 
(log;y Re — 2) 


During a water-takeoff run for an amphibious aircraft without a hydrofoil, this viscous 
resistance term dominates the drag in water during the planing stage. In the model 
considered, when the hull is unported from the water with the help of the hydrofoil, viz. 
the buoyant force is zero, this term is also zero. 

The kinematic equations are modeled via a time-stepping approach, in which the 
speed u, time t, and displacement x variables are discretized, with subscripts 0 <i < N; 
for N; timesteps. The dot notation is used to denote time derivatives. This approach is 
similar to the water-takeoff model for amphibious aircraft used in the work of Qiu and 
Song [8]; they did not, however, consider the implementation of a hydrofoil in their study. 
The following equations are solved for the speed and displacement by integration of the 
equations of motion for some initial condition (to, Xo, Ug, Uo). The output consists of 
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variables including the water-takeoff distance xy, the time taken tw, and the other relevant 
parameters for analysis and optimization. 


Uj = Uj_y + Uj (t; — ti-1), (14a) 
1 
Xj = Xj-4+ 5 (ui Uj) (E— 44). (14b) 


In the water-takeoff regime, the wetted area of the hydrofoil underwater with dihedral 
angle 0, ¢ (negative, hence anhedral) reduces as the waterline height of the aircraft decreases 
due to the lift generated by the hydrofoil. The lift generated by the submerged area of the 
hydrofoil with cavitation effects as an initial approximation can be modeled as: 


1 
le 5 PWH CL Sf (I), (15) 


where the submerged area Sj... (4), now a function of the height , is assumed to be fully 
wetted underwater, and the lift coefficient of the hydrofoil takes the losses due to cavitation 
effects into account. For an initial approximation assuming the waterline interface as a flat 
surface, this underwater wetted area can be determined geometrically as: 


2 
Shf, h 

S h) & <1 by + —— |, —90° < by¢ < 0°. 16 

hf wee (1) bP, ( hf a) hf (16) 





To obtain the lift and drag forces over the hydrofoil, either experiments or simulations 
including at least cavitation modeling need to be performed to obtain a reasonable estimate 
at this stage of the design process. It is computationally very expensive to determine the 
hydrodynamic forces via high-fidelity CFD at every time-discretized point in a takeoff 
analysis as observed by Seth and Liem [44], so a surrogate model, or an approximation 
model, is generated for use in the water-takeoff analysis for the optimization problems to 
reduce the computational burden. The surrogate models will compute the non-dimensional 
coefficients of the hydrofoil as functions of speed and angle of attack. 

In this work, sample-based surrogate models are selected for their simplicity and 
non-intrusive nature. To construct the surrogate models, we first need to generate samples 
by running the high-fidelity analyses. A design of experiments is performed to generate 
the design space by selecting N in the speed range ug < u < U7o, with initial speed up, to 
generate a set {u} and Ny samples in the angle of attack range amin < & < Qmax to generate 
a set {a}. The Cartesian product of these sets defines the design space to generate CFD 
training and testing data for the surrogate, as will be discussed in Section 4.2.2. Similar 
arguments can be applied to airfoils, so surrogate models are also considered for them in 
the framework, although they may not be necessary with the use of low-fidelity solvers, 
which is to be further discussed in Section 3.4.2. The water-takeoff analysis with surrogate 
models is presented as an XDSM in Figure 4. 

Note that this takeoff procedure alone does not guarantee the aircraft will satisfy the 
required stability and trim characteristics, which are considered in the following section to 
further elaborate the considerations for the optimization problem presented in Table 2. 
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Figure 4. XDSM diagram for the water-takeoff simulator with surrogate models. 


3.3. Stability Considerations and Analysis 


In this section, we first discuss the physical considerations in evaluating the longitudi- 
nal stability requirements of amphibious aircraft, followed by the details of trim analysis. 
As we assume a no-wave condition, the lateral stability conditions caused by irregular 
wave patterns and effects of sponsons during water-takeoff are not considered at this 
preliminary design stage. 


3.3.1. Physical Considerations 


Consider Cartesian coordinates indicating the position of the aerodynamic center of 
the hydrofoil (xj,,Znf) where the origin is set at the aircraft CG, assumed to be located 
at 33% MAC of the wing aft of the leading edge for the following analysis. The large lift 
and drag forces produced by the hydrofoil will generate significant moments. Looking 
at an amphibious aircraft with its nose facing left, the hydrofoil will be located below 
the CG (thus, z;r¢ < 0), and its drag will resultantly contribute a nose-down pitching 
moment, which must be corrected by the elevator to maintain the required hull trim angle. 
To minimize elevator effort, it is best to counter this anti-clockwise moment by placing 
the hydrofoil ahead of the CG (thus, x;¢ < 0), so the hydrofoil is able to provide lift that 
also provides an opposing pitch-up moment to compensate for the moment caused by its 
drag. This motivates a positioning procedure, such that the minimum takeoff distance is 
achieved during the water-takeoff run while maximizing moments generated by lift and 
minimizing moments generated by drag. A complication of this procedure is that the lift 
and drag of the hydrofoil suffer from variations in speed and angles of attack dissimilar 
to their airfoil counterparts due to cavitation and ventilation effects and are, thus, not 
predictable via tools meant for airfoil analysis; hence, the considerations using high-fidelity 
CFD and surrogate models introduced in Section 3.2. The effects of ventilation, however, 
are not explicitly modeled in this work, for reasons described in Section 2.2. 


3.3.2. Trim Analysis 


The non-dimensionalized moment equation is derived in the body-axis system of the 
aircraft following the process of Raymer [45] with modifications for the hydrofoil. The 
drag of the hydrofoil is not considered as negligible in this derivation unlike as usually 
considered for the tail, as the density of water is approximately three orders of magnitude 
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ereater than that of air. The following expression is solved for the trim condition of the 
moment coefficient about the center of gravity Cy, = 0, to determine the required C_,, to 
correct the moments generated by the hydrofoil and the wing. 


Chg = Cm, + CMeng Cr. (Xcg — Xac,) — MnVnCr, + Ynf Ving Cop — Zag Coy, |, (17) 


where C Vee is the moment coefficient corresponding to the engine, %c¢g = Xcg/Cw is the 
non-dimensionalized location of the center of gravity and Xqc,, = Xac, /Cw is the location of 
the aerodynamic center of the wing with respect to the nose of the aircraft as the origin. 
The dynamic pressure ratios 1, and 1, are defined as: 


1 y) t 2 
ZP AU), tnt = 2PW" nF 
i / hf = 4 
sP Aus, sP Aus, 





qh = (18) 


The longitudinal stability coefficients for a component c corresponding to a lifting 
surface are defined as: 


Vo =a-z Le re (19) 


with moment arms |,¢ = Xcg — Xn, Zhf = Zcg — Zn¢, Which will be determined via the 
optimization procedure for the hydrofoil location presented in Section 3.1.4. 


3.4. Design Framework and Solvers 


In this section, we summarize the entire framework of the hydrofoil design process 
described in the previous sections, and present the relevant computational solvers for 
analyses. 


3.4.1. Overall Design Process 


The entire design process is depicted in the XDSM presented in Figure 5, with the 
descriptions of diagram components are as described in Section 3.1.3. The water-takeoff 
simulator is first executed without any hydrofoil to obtain the initial aircraft design’s water- 
takeoff performance. The hydrofoil design parameters, determined from the preliminary 
sizing, are then used for the generation of its surrogate model using high-fidelity CFD, and 
similarly for the wing if required, for the design space of takeoff speeds and operational 
angles of attack. These surrogate models are used in the minimization of the water-takeoff 
distance to determine the optimal span and incidence angle of the hydrofoil. These optimal 
values are then provided to the optimization loop to minimize the required stabilizer force 
to achieve trim during the water-takeoff run by adjusting the position of the hydrofoil 
for the optimum. The relevant data, such as the optimized takeoff distance, span length, 
incidence angle, minimum stabilizer force, hydrofoil location, and the data obtained from 
the water-takeoff procedure, are provided to the designer to assess the performance and 
feasibility of the design. 


3.4.2. Computational Fluid Dynamics Solvers 


The choice of the computational method of CFD for the determination of the coeffi- 
cients of the airfoil and hydrofoil is extremely important, depending on the appropriate 
compromise between accuracy and speed at the required stage of the design process. 
Reasonably accurate low-fidelity solvers with good computational efficiency are openly 
available for airfoil analyses, but the same are not openly available for hydrofoil analyses 
with cavitation models, so open-source high-fidelity solvers are considered instead. In this 
work, low-fidelity solvers are deemed sufficient to evaluate the aerodynamic performance 
of aircraft wing within the operational regime considered in this study, while high-fidelity 
solvers are required to evaluate the hydrodynamic performance of hydrofoils, as described 
briefly below. 
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Figure 5. XDSM diagram for the hydrofoil design process. 


3.4.2.1. Low-Fidelity Solvers 


To reduce computational burden, XFOIL [41], a well-established viscous-inviscid 
design and analysis tool for airfoils at low Reynolds numbers, can be used to compute the 
aerodynamic coefficients for the wing by computing the lift coefficients of the airfoil Cy, 
with flap AC; and ground effect AC), corrections in the takeoff simulation itself: 


Cre = Cr, a ACL, zi ACy,. (20) 


Note that XFOIL may not be able to predict the coefficients for flows that occur at 
high angles of attack beyond the stall angle, in which case high-fidelity CFD must be used; 
however, angles greater than the stall angle are not usually reached until the aircraft rotates 
after reaching the takeoff speed. A surrogate model via high-fidelity CFD can also be 
generated for the airfoil to account for ground effects more accurately during takeoff. 


3.4.2.2. High-Fidelity Solvers 


The unsteady analyses with cavitation modeling over the hydrofoil are performed by 
solving the URANS equations via the interPhaseChangeFoam module of the OpenFOAM 
framework [46], following the approach of Vernengo et al. [29]. One difference from this 
analysis is that the Reynolds numbers for these cases are in the turbulent regime; hence, the 
Spalart-Allmaras turbulence model is used based on the approach of Garg et al. [28]. The 
cavitation model used in this analysis is the Schnerr-Sauer model [47], and the equations 
are solved using the segregated PIMPLE approach. 


4. Case Study Description and Results—10-Seater Amphibious Aircraft 


A preliminary design of a 10-seater aircraft based on a DHC-6 Twin Otter is considered 
for the case study. A planing hull from the National Advisory Committee for Aeronautics 
(NACA) Technical Note TN-2481 report [49] is adopted for this design, where we derive a 
curve-fit from the available experimental data for the water-takeoff simulation purpose, as 
shown in Figure 6. The aircraft weight is approximately W = 5620 x 9.81 N, its wing area 
is approximately S,, = 39.2 m* with an incidence angle of a» = 0°. The takeoff speed is 
selected to be 44 m/s, similar to the Twin Otter’s. The selection of the airfoil for the wing is 
the NACA 63(4)-12 profile. A surrogate model for this airfoil using RANS without ground 
effects has been generated by Seth and Liem [1], which is used here. 
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Figure 6. NACA TN-2481 planing hull data. 


A quadratic thrust model presented by Gudmundsson [50] is used based on the 
specifications of a Pratt and Whitney Canada PT6A-34 turboprop engine via the follow- 
ing expression with T as the thrust of the engine as a function of speed u, as shown in 
Equation (21). Figure 7 shows the thrust model only within the takeoff speed range of 
interest; note that the bounds of the domain make it look more linear than quadratic. 
The engine thrust is scaled for the first 10 s of the takeoff according to the expression 
in Equation (22) as a heuristic pilot model for the increase in thrust, similar to the case 
considered in Gudmundsson [50]. 


Totatic — 21 T, — 2T stati 
T(u) = ( static : tom) 42 Be (= Umax sate u+ ae (21) 
Umax Umax 
0.25 +0.75(t/10) £< 10, 
Ce (22) 
0 t > 10. 
13,000 
12,000 
11,000 
10,000 
9,000 
0 10 20 30 AO 


Velocity, m/s 


Figure 7. Thrust model as a function of speed, based on the specifications of a Pratt and Whitney Canada PT6A-34 
turboprop engine. 
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4.1. Sizing and Profile Selection 


The elevator effectiveness speed is considered to be uj, ~ 15 m/s. The wing lift at 
this speed with the corresponding trim angle at aim = 6° is estimated to be C;,, = 1.5 
using XFOIL and assuming AC; + AC;, = 0.3. This does not account for interference 
effects between the wing and the fuselage. The hydrofoil profile selected for this case is the 
Waid-Lindberg hydrofoil profile [40] in a strut-based configuration. The cavitation number 
at this speed is Ca © 0.87, which is in the cavitating regime for the hydrofoil. The C,,, 
in this regime according to the experimental data is approximately 0.9, and the hydrofoil 
size can, hence, be determined by Equation (8), giving Sj¢,,, © 0.455 m*. Considering a 
retractable, non-foldable hydrofoil system, an aspect ratio of 6 is selected initially, giving a 
chord length ¢;,¢ ~ 0.275 m and initial span value bj, + 1.65 m. The relevant parameters 
for the study are tabulated in Table 3. 


Table 3. Aircraft and hydrofoil parameters. 


Parameter Value 

Wing Area, m* 392 
Wing Aspect Ratio 10.0 
Wing Span, m 19.8 
Wing Chord, m 1.98 
Wing Stall Speed, m/s 36.8 
Wing Incidence Angle, ° 0.0 
Wing Dihedral Angle, ° 3.0 
Hull Beam Width, m 1.71 
Takeoff Speed, m/s 44.0 

Hydrofoil Area, m? 0.454 
Hydrofoil Aspect Ratio 6.0 
Hydrofoil Span, m 1.65 

Hydrofoil Chord, m 0.275 
Hydrofoil Mass, kg 60.0 


4.2. CFD and Surrogate Model Generation 


This section discusses the use of CFD and surrogate models in the case study. The 
erid generation, CFD analyses, and surrogate model validations are described. 


4.2.1. Grid Generation, Convergence Study, and Validation 


The meshes are generated using a hyperbolic grid generator via the Python mod- 
ule pyHyp (https://github.com/mdolab/pyhyp) developed by the Multdisciplinary De- 
sign Optimization Lab (MDOLab) at University of Michigan based on the theory from 
Luke et al. [51]. This generates structured hexahedral meshes which keep a low cell count 
with high quality, which are then converted into the unstructured OpenFOAM format for 
this study. The domain extends to approximately 30 chord lengths from the hydrofoil. The 
sharp leading edge of the hydrofoil is modified to create a small, blunt leading edge to 
ensure high quality grid generation, and a large number of elements is concentrated near 
the blunt trailing edge to resolve unsteady behaviour more accurately. At this prelimi- 
nary design stage, we consider non-accelerating flows using a constant speed boundary 
condition. 

The meshes and solver are validated using the experimental results from Waid and 
Lindberg [40] at Re = 7.8 x 10°, Ca = 0.293, u =9.144 m/s, ene = 6°, and three grids (Lo 12) 
are analyzed via Richardson extrapolation of the time-averaged drag coefficient Cp, yr using 


a grid refinement ratio of ,/ Nj_1/N; = 2, where N; is the number of cells of the coarsened 
erid with respect to the finer grid N;_;. A wall distance of 2 x 10~* m is selected to match 
y* > 30 with the use of wall functions using the Spalart-Allmaras turbulence model to 
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avoid extremely small timesteps, which would drastically increase the computation time. 
The appropriate chord length cj and non-dimensional time T are selected such that the 
cavitating flow is developed over the entire chord length of the hydrofoil to determine the 
required simulation time t, given by Equation (23): 


t 
ee, (23) 
Chf 


The L; grid and the leading edge modification mentioned previously are shown in 
Figure 8. The extrapolation is shown in Figure 9. The results are provided in Table 4, with 
a depiction of the cavity over the hydrofoil shown in Figure 10. Computing the relevant 
ratio of the grid convergence indices (GCI) for Cp, ,; as the quantity of interest gives 
approximately 0.944 ~ 1 within the asymptotic range and approximate solver convergence 
p = 1.34. The error from the L, grid is considered as acceptable for a preliminary design 
study with lower computational costs of analyses. The analyses at this Reynolds number 
are within the laminar regime, and the use of turbulence models may have contributed 
to the error. Further analyses on this grid are performed at higher Reynolds numbers 
in the water-takeoff regime for amphibious aircraft, for which experimental data are 
not available; hence, this design point was used to validate the convergence study with 
available experimental data. 
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Figure 8. L; grid resolution and blunt leading edge with velocity contours (in m/s). 
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Figure 9. Grid convergence study via Richardson extrapolation. 
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Water Volume Fraction 


Figure 10. Supercavity over Waid-Lindberg profile. 


Table 4. Grid convergence study and validation. 


Grid Grid Size Cop, ; Error Cr ; Error 
Experiment - 0.0421 0% 0.4199 0% 
Lo 3950 0.0522 23.99% 0.4667 11.14% 
Ly 15,958 0.0458 8.79% 0.4556 8.5% 
Lo 63, 832 0.0432 2.61% 0.4480 6.69% 


4.2.2. Surrogate Models 


The surrogate models for time-averaged Cy, Cp, and Cy of the supercavitating hydro- 
foil are presented in Figure 11. Note that the profiles shown here significantly differ from 
ones generated by airfoils, as the cavitation effects are accounted for [1]. The surrogate 
models are generated using Kriging models via the SMT: Surrogate Modeling Toolbox [52] 
for the Waid-Lindberg hydrofoil profile using 25 sample points over 0 < ayj¢ < 15 and 
2 <u < 45. The 25 samples are selected from a data-set of uniformly distributed 100 points 
generated via the high-fidelity computational model mentioned in Section 3.4.2.2, i.e., the 
unsteady analyses with cavitation modeling. The timestep set for the CFD analyses was 
2 x 10~° s, and the iterations indicated a maximum Courant-Friedrich-Lewys (CFL) num- 
ber of 2 per timestep. The remaining 75 samples, denoted by x;,, are used for validation 
purposes, in which we compute the normalized root-mean square deviation (NRMSD), as 
shown in Equation (24). The errors are normalized by the ranges of the corresponding non- 
dimensional coefficients. The computed NRMSD values for C,, Cp, and Cy are tabulated 
in Table 5. These errors are deemed acceptable in approximating the coefficients. 


R = Xmax — Xmin, Ci = Xi — Xgin XE {Cp, Cr, Cu}. (24) 
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The (L/D)n¢ = (CL/Cp)ng ratios, computed pointwise over the design space, are 
presented as a surface in Figure 12. The drastic reduction of (L/D)j,¢ at higher speeds, 
which is caused by cavitation, is clearly visible in this profile. 


Table 5. Error analysis of surrogate model for Waid-Lindberg supercavitating hydrofoil. 


Cp Cr Cm 
NRMSD, % 2.8212 3.0168 3.8037 
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Figure 11. Surrogate contour plots for the coefficients of the Waid-Lindberg supercavitating hydrofoil. These profiles 
include the cavitation effects. 
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Figure 12. Response surface for the Waid-Lindberg supercavitating hydrofoil’s (L/D) r surrogate. 
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4.3, Optimization Studies 


The optimizer used for the following studies is the Sequential Least-Squares Quadratic 
Programming algorithm provided via the SciPy package in Python [48]. This optimizer is 
suitable for optimization problems with bounds and constraints. The optimizers are set to 
converge at a tolerance of 107. 


4.3.1. Water-Takeoff Distance Minimization 


A water-takeoff distance minimization study is performed to find the optimum span 
and incidence angle of the hydrofoil, following the procedure described in Section 3.1.3. 
As a retractable, non-foldable configuration is considered, the span is bounded from above 
by the hull beam width for the hydrofoil system to fit into the hull upon retraction. 

The constrained design space with the baseline and optimum points are shown 
in Figure 13, in which dark dots with larger radii indicate design points in the feasible 
region at which the lift constraint is satisfied, and the colors corresponding to the colorbar 
depict the takeoff values for all points. This diagram also depicts the sensitivities of the 
objective function to the design variables. The number of iterations required to obtain 
convergence was 16, with 181 function evaluations and 16 gradient evaluations. The 
water-takeoff distance exhibits nonlinear behavior with respect to both variables, and the 
minimization problem appears to be multimodal and sensitive to initial guesses. Observing 
slices of fixed span lengths with varying incidence angles indicates that the takeoff distance 
sharply increases at specific incidence angle values. As the span length increases, the 
corresponding incidence angle decreases, and eventually extends to multiple incidence 
angles as can be seen in the region 2.2 m < by, _ 4.0m, 0° < a, fs 4.0°. This physically 
corresponds to the regions in which the hydrofoil generates large drag and insufficient 
lift to reduce the hull resistance effectively. At larger angles, such as the region between 
6-14° for span lengths between 1m < by¢ < 6m, the hydrofoil generates larger lift, hence 
reducing the takeoff distance values; however, this region violates the lift constraint and 
the optimizer chooses the lower takeoff distance values satisfying the constraint at the 
optimum shown in the diagram. 
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Figure 13. The design space for water-takeoff distance minimization to find the span and incidence angle of the hydrofoil. 


The darker region with larger radii depicts the area satisfying the lift constraint. 
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If an incidence angle of #;, = 5° is selected as an initial guess with the initial span 
value, then the optimizer tends to miss the minimum at lower angles, as it tends towards 
higher angles. By choosing the right starting point, on the other hand, the global optimum 
within the bounds can be found more effectively. In particular, the initial guess is selected 
through an initialization process that aims to maximize the average (L/D), f- Phe maxi- 
mum average (L/D)j¢ ratio over the takeoff speed range of the hydrofoil is evaluated as 
an initial guess to feed to the water-takeoff distance minimizer. This is done by evaluating 
the averages of the (L/D), diagram in Figure 12 over the speed range, providing the 
average (L/D)n, ratio as a function of the angle of attack of the hydrofoil. This procedure 
results in the graph shown in Figure 14. The maximum point obtained in this initialization 
procedure is then used to determine the incidence angle’s initial value in the water-takeoff 
distance minimization procedure. A similar profile was obtained from a study on the 
YS-920 hydrofoil, designed for subcavitating conditions [1]. However, it suffered from 
lower (L/D)nf ratios and higher water-takeoff distances, indicating profiles designed for 
subcavitating regimes were not suitable for amphibious aircraft and further emphasized 
the need for using supercavitating hydrofoils in this particular application. 


@ Baseline 


0 2 4 6 8 10 12 14 


Angle of Attack ap, r, ° 


Figure 14. Average (L/D)},¢ Ratios for the Waid-Lindberg supercavitating hydrofoil. 


4.3.2. Stabilizer Force Minimization 


The optimum span and incidence angle are then provided to the stabilizer force 
minimization problem, to find the optimal hydrofoil position. The optimization procedure 
follows the description presented in Section 3.1.4. The initial values selected for the 
hydrofoil location are based on the criteria presented in [3] that the longitudinal position 
of the hydrofoil be located at approximately 0.5¢,, forward of the aircraft’s CG. Its design 
space is depicted in Figure 15, which also depicts the sensitivities of the objective function to 
the design variables. It indicates that the stabilizer force is more sensitive to the horizontal 
location of the hydrofoil than its vertical location. It also indicates that the upper bound for 
Zhfr shown as the horizontal red line in the figure, will be the solution, and will calculate 
the appropriate x;,¢ to minimize the stabilizer force, which is physically consistent with 
the considerations presented in Section 3.3.1. The number of iterations required to obtain 
convergence was 18, with 98 function evaluations and 18 gradient evaluations. 
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Maximum Stabilizer Force, N 
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Figure 15. Hydrofoil position design space for the stabilizer force minimization. The red line indicates the upper bound of 


the hydrofoil’s vertical position. 


The results of both optimization procedures are summarized in Table 6. Detailed 
discussions on the performance comparison between the baseline and optimized designs 
are presented in the next section. 


Table 6. Span length, incidence angle, and position optimization results. 


Optimization Problem Parameter Optimal Value 
Distance, xw 359 m 
Water-Takeoff Distance Incidence angle, ay ¢, 1.486° 
Span, Dig 0.914 m 
Maximum Stabilizer Force H nse Lue — 
aan (—0.5959 m, —2.48 m) 
(xp fr hf ) 


4.4, Performance Evaluations, Comparisons, and Discussion 


In this section, we evaluate the performance of hydrofoils in the context of amphibious 
aircraft application. The water-takeoff analysis procedure described in Section 3.2, yields 
detailed information of force and load variations during takeoff, which provides insights 
into the effectiveness of adding hydrofoils. 

First, we compare the water-takeoff performance of the aircraft with and without hy- 
drofoils; note that the optimized hydrofoil configuration is considered in this comparison. 
The comparison of its effects on the hull resistance is shown in Figure 16. It indicates minor 
reductions in the wave resistance of the hull, and substantial reductions in the viscous resis- 
tance. These observations physically correspond to the hydrofoil reducing the submerged 
volume and wetted area during the water-takeoff by the hull being unported from water, as 
can be seen in the water-takeoff results shown in Figure 17, satisfying the purposes of the 
design. The viscous resistance of the configuration without hydrofoils indicates the hull is 
still submerged to some extent as it reaches the takeoff speed, whereas hydrofoils can help 
completely unport the hull before this speed. The optimization procedure of the incidence 
angle and span ensures that the lift constraint in Table 1 is satisfied; the hydrofoil’s lift 
reduces as the wing’s lift increases while the aircraft’s speed is below the required takeoff 
speed, as observed in Figure 18. 
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Figure 16. Hull resistance comparisons with and without the hydrofoil. 
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Figure 17. Variations of hydrodynamic and propulsive forces during water-takeoff process. 


The results of the aircraft configuration with the baseline hydrofoil design (i.e., when 
using the initial values of the optimization design variables) are compared to those with- 
out and with the optimized hydrofoil design in Figure 19 to determine the performance 
differences and the relevant physics. A drawing depicting the differences between the 
baseline and optimized designs is presented in Figure 20. The comparisons between the 
three configurations in terms of water-takeoff distance and maximum stabilizer force are 
tabulated in Table 7. The notable reductions of both objective functions in the optimized 
design, compared to the baseline design, are evident in these results, which confirm the 
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effectiveness of the optimization procedure. More importantly, the results show that the 
optimized hydrofoils help reduce the water-takeoff distance by approximately 3.5% and 
the maximum stabilizer force by almost half as compared to those of the no-hydrofoil con- 
figuration, demonstrating the benefits of adding hydrofoils to amphibious aircraft. These 
results also emphasize the importance of performing optimizations, as adding hydrofoils 
based on physical reasoning alone does not guarantee performance improvement. The 
water-takeoff distance achieved by the aircraft configuration with the baseline hydrofoil, 
for instance, is longer than that of the no-hydrofoil configuration, which would defeat the 
purpose of adding hydrofoils. 
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Figure 18. Lift and drag forces of components during water-takeoff process. 
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Figure 19. Force comparisons of the amphibious aircraft with and without the hydrofoil, considering both the baseline and 


optimized hydrofoil configurations. 


The baseline design achieves the aim of reducing the buoyant load, shown in Figure 19, 
before the no-hydrofoil configuration, but it generates excessive moments that may not be 
corrected by the elevator, based on the required stabilizer force also seen in Figure 19. The 
relatively lower absolute values of the stabilizer force required for the optimized design 
as compared to the baseline design indicate the optimization procedure for the hydrofoil 
location is appropriate. Physically, the optimized design exhibits growing sinusoidal 
oscillations of the required stabilizer force, indicating a porpoising tendency for the aircraft 
during the water-takeoff if the trim corrections are not performed. 
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Figure 20. Aircraft drawing with baseline and optimized hydrofoil configurations (in meters). 


Table 7. Water-takeoff performance results and comparisons between the different configurations. 


No Hydrofoil Baseline Optimized 
Water- Takeoff Distance, m 372 461 359 
Maximum Stabilizer Force, N 604 4743 435 


5. Conclusions 


In the present study, the key considerations for the preliminary design of a hydrofoil 
for amphibious aircraft were investigated. A sizing procedure was formulated with recom- 
mendations for profile selection and different configurations. In particular, we observed 
that a supercavitating hydrofoil profile was required for amphibious aircraft application. 
An optimization framework was formulated to minimize the water-takeoff distance of the 
aircraft by finding the optimal span length and incidence angle. Its detrimental effects on 
the longitudinal stability of the aircraft with respect to the required horizontal stabilizer 
force to maintain the trim angle were also minimized, by finding the optimal position of 
the hydrofoil with respect to the aircraft’s center of gravity. The effectiveness of the design 
was analyzed by a water-takeoff physics simulator using surrogate models of the hydrofoil. 
These models were generated via CFD sample data to determine the performance of a 
selected hydrofoil over a range of speeds and angles of attack within the water-takeoff 
regime below the takeoff speed. 

To evaluate and assess the performance and effectiveness of adding hydrofoils to am- 
phibious aircraft, the developed preliminary design framework was applied in a 10-seater 
aircraft as a case study. Results from the water-takeoff simulator showed the variations 
of forces and loads during the water-takeoff run, providing insights into the physics and 
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behavior of the aircraft during the process. The aircraft’s porpoising tendency, for instance, 
was evident by observing the stabilizer variations, which highlighted the importance of 
trim corrections. These results were unique to amphibious aircraft, i.e., the hydrofoil 
performance was different from that of marine applications, and the takeoff performance 
was different from that of conventional ground-based aircraft. 

The results of the water-takeoff procedure indicated that the hydrofoil performed 
the required purposes of reducing hull resistance and buoyant load. These reductions 
demonstrated the hydrofoil’s effectiveness in improving the water-takeoff performance and 
efficiency of amphibious aircraft, as evident in the reduction of the water-takeoff distance 
observed in the case study. Our results showed, however, that adding hydrofoils alone, 
without being properly optimized, could not guarantee a performance improvement. This 
observation further emphasized the need to develop a comprehensive design framework 
for amphibious aircraft and highlighted the benefits of performing optimization in design 
studies. The optimization results indicated that the water-takeoff distance minimization 
with respect to the span and incidence angle as design variables was non-differentiable 
and multimodal with the required constraints, depending on the hydrofoil’s profile. The 
problem of minimizing the stabilizer force was relatively less complex, but still an important 
part of the design procedure. 

To the best of our knowledge, ours is the first study that systematically evaluates 
hydrofoil performance in the context of amphibious aircraft operation, particularly during 
the water-takeoff process. The modular design of the framework gives flexibility to replace 
the current models, be it aerodynamic, hydrodynamic, or propulsive, with higher-fidelity 
models in future development. The computation with the surrogate models allows for 
rapid analysis with variations of parameters, taking less than two minutes to complete 
for a case, which offers a desirable computational efficiency. This computational frame- 
work will reduce reliance on expensive experiments which have limited the technological 
advancement and progress in amphibious aircraft development. More importantly, the 
developed framework will enable performing detailed amphibious aircraft design opti- 
mization and tradeoff studies with more computational rigor. We will be able to perform 
multidisciplinary design optimization that simultaneously considers the aerodynamic and 
hydrodynamic performance (e.g., to optimize hull and hydrofoil shapes) to obtain a truly 
optimal design, thereby advancing the amphibious aircraft development. 
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Abbreviations 


The following symbols used in the manuscript are provided here for reference: 


a" ™N 


XZ) 


——™~ 
VY 


Poo, Patm, PV 


2D coordinates with respect to center of gravity as origin 
Time 

Non-dimensional time 

Time-averaged quantity 

Change in quantity 

Time derivative of quantity 

Acceleration due to gravity on Earth 

Reference height for hydrostatic pressure head 

Speed in horizontal direction 

Freestream, atmospheric and vapor pressure 


re Density of fluid 
i Dynamic pressure ratio 
u Dynamic viscosity of fluid 


Re, We, Fr, Ca 
Ow, One Ong 


Reynolds, Weber, Froude, and cavitation number 
Quantity corresponding to wing, horizontal tail, and hydrofoil 


L, D, M Lift, drag and moment 
Ci Cpe Cu Lift, drag, and moment coefficients 
Cp, Cr Pressure and skin-friction resistance coefficients 
Ca, Cr, Cy Load, resistance and speed coefficients 
Wi Dey Gin Angle of attack, incidence, and trim angle 
L/D Lift-to-drag ratio 
Ve, Le Longitudinal stability coefficients of component c 
Leeda Horizontal and vertical moment arms 
bic Span and chord lengths 
S Reference area 
L Reference length 
AR Aspect ratio 
WwW Aircraft weight 
B Hull beam width 
d Dihedral angle 
C Mean aerodynamic chord length 
(X, Z) Coordinates non-dimensionalized with respect to mean aerodynamic chord 
T Thrust 
QT Thrust angle 
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Abstract: The design of aircraft and engine components hinges on the use of computer aided 
design (CAD) models and the subsequent geometry-based analyses for evaluation of the quality 
of a concept. However, the generation (and variation) of CAD models to include radical or novel 
design solutions is a resource intense modelling effort. While approaches to automate the generation 
and variation of CAD models exist, they neglect the capture and representation of the product’s 
design rationale—what the product is supposed to do. The design space exploration approach 
Function and Geometry Exploration (FGE) aims to support the exploration of more functionally 
and geometrically different product concepts under consideration of not only geometrical, but also 
teleological aspects. The FGE approach has been presented and verified in a previous presentation. 
However, in order to contribute to engineering design practice, a design method needs to be validated 
through application in industrial practice. Hence, this publication reports from a study where the FGE 
approach has been applied by a design team of a Swedish aerospace manufacturers in a conceptual 
product development project. Conceptually different alternatives were identified in order to meet 
the expected functionality of a guide vane (GV). The FGE was introduced and applied in a series 
of workshops. Data was collected through participatory observation in the design teams by the 
researchers, as well as interviews and questionnaires. The results reveal the potential of the FGE 
approach as a design support to: (1) Represent and capture the design rationale and the design space; 
(2) capture, integrate and model novel solutions; and (3) provide support for the embodiment of 
novel concepts that would otherwise remain unexplored. In conclusion, the FGE method supports 
designers to articulate and link the design rationale, including functional requirements and alternative 
solutions, to geometrical features of the product concepts. The method supports the exploration of 
alternative solutions as well as functions. However, scalability and robustness of the generated CAD 
models remain subject to further research. 


Keywords: product development process; design space exploration; function modelling; 
CAD; aerospace 


1. Introduction 


The development of new components for aircraft engines often relies on improvement and 
refinement of an existing “legacy design”. To be able to meet the targets set by Flightpath 2050 [1], 
as well as contemporary developments in air travel [2], manufacturers need to explore and develop 
conceptually new designs. Such novel concepts are required on an aircraft and engine system level, 
but equally on a component level [3,4]. 
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Examples for such novel product concepts are the open rotor design [5] or electrical and hybrid 
propulsion designs [6]. Such concepts bring new arrangements of parts, and require new functions to 
be met by the engine components. 

Such radically new concepts—as well as the components they are composed of—need to be 
developed, designed and evaluated first. In the case of developing physical products, such as 
aeroengine components, this requires their embodiment in a sufficient level of detail for engineering 
analysis. In most cases, this means that a geometry model in the form of computer aided design (CAD) 
of each concept needs to be created. To be able to assess all relevant product properties, analyses on 
solid geometry models may be required [7]. As a result, this publication is focused on the management 
and generation of solid CAD models. 

To be able to develop novel products outside the already known design space, developers 
need reliable methods and tools that support them in the generation of concepts and models as 
well as their analysis and evaluation [8-10]. Currently available—and applied—methods such as 
parameterisation [11] or knowledge-based engineering (KBE) [12] cover only a dimensional variation 
of existing concepts or require an expensive setup of a master model. Therefore, these geometry based 
approaches only allow for the exploration of a subset of the available design space [8,13]. To enable 
conceptual design, which explores alternative functions and solutions, a design space exploration (DSE) 
method needs to take into account a representation of the product’s function [14,15]. 

The DSE method function and geometry exploration (FGE) has been developed to fulfil this 
premise: through a connection of the functional and geometrical domain, a wider and more efficient 
exploration is to be enabled [16]. However, while the FGE approach itself has been presented and 
verified in [16], its effect on improving DSE in practice has not been validated. This publication reports 
from a study to validate the approach through application in an industrial setting. Validation refers in 
this case to test whether it is useful, that means it has a positive impact on the product development 
process in an applied context [17,18]. 

The validation is performed through a qualitative study about the application of the FGE approach 
in a development project for a static turbine structure, a fan frame guide vane (GV). The study has 
been conducted in tight collaboration with an aerospace manufacturer. The study investigates into 
whether FGE does improve DSE in that it enables the capture, representation and embodiment of 
novel design concepts. 

An over view over the underlying theories is given in Section 2, and the FGE approach itself 
is shortly explained in Section 2.1. How the study has been setup and how the approach has been 
evaluated is explained in Section 3. The application of FGE in the industrial setting is presented in 
Section 4. The practitioners’ evaluation and reflections on the application of the approach are presented 
in Section 4.2. Section 5 discusses these results in relation to the problematic described above as well 
as compared to similar approaches from literature. Lastly, Section 6 summarises the contribution of 
this paper and proposes directions for further development of the approach. 


2. Frame of Reference 


In the conceptual phases of product development, multiple alternative concepts need to be 
evaluated. The sum of all possible concepts which fulfil the respective requirements is called the design 
space [19]. The more, and more different, concepts can be investigated—corresponding to a larger area 
of the design space being covered—the higher the probability that a good, or even the best possible, 
concept has been selected [8]. 

To achieve this, the design space needs to be defined and the concepts residing in it need to be 
devised and evaluated against each other [9]. For the evaluation against each other, that is assessing 
whether, and how well, they fulfil requirements and functions, the concepts’ need to be simulated in 
terms of different behaviour. For most of such simulations, a geometrical product representation in the 
form of a CAD model is needed [10]. 
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CAD models are a representation of a product’s geometry and topology. Geometry refers to the 
shape of the product, whereas topology refers to the relations of the individual shapes towards each 
other [20]. Most commercial CAD software operate as feature-based parametric CAD, where geometric 
and topological information is contained in features, which are set into relation in a procedurally 
operated feature-tree. If set-up and maintained correctly, the alteration and exchange of such features 
should allow for easy adjustment of the CAD model [21]. Furthermore does parameterisation allow 
for a variation of geometry model’s dimensions [22]. 

While parameterised, feature-based CAD models are theoretically easily altered, in practice the 
alteration of is difficult. A lack of consistent modelling practices, the complexity of relations in a master 
model and ambiguous feature definitions make CAD models difficult to edit and lead to failures in 
their regeneration [20,21,23]. Furthermore is parameterisation considered “a big challenge”, and even 
in the best case “offers little flexibility”, since the parameterisation of the master model already defines 
all possible concepts and does not allow for the introduction of novel solutions [24]. 

To avoid such failures, CAD models are build following design guidelines. While the composition 
of a CAD model’s structure represents the model’s design intent, this does rarely coincide with the 
product’s design rationale (DR). Hence are changes to the product, which are based on the introduction 
of new functions or new solutions for identified functions, difficult, since the CAD model does not 
contain any of this teleological information [14,15]. 

Function models on the other hand can represent a product’s DR, together with architectural or 
systemic product properties. As such, they provide a good basis for teleologically motivated DSE: 
through the representation of functions and solutions, developers can precisely target their design 
effort [25]. Furthermore can function models be used for systemic analysis [26], automated DSE 
and concept generation [27,28] or product platforms [25,29]. However, while in some cases enabling 
“physics-based reasoning” [30], function models do not provide a basis for most engineering analysis, 
where product geometry models are needed. Connecting the two modelling domains of function and 
geometry has been a long standing challenge in engineering design research (EDR) [14,15,31]. 

enhanced function-means (EF-M) modelling is a function modelling approach which represents 
the design space through alternating functional requirement (FR) and design solution (DS) objects [32]. 
A DS can be subject to individual constraints, which place limits on design parameters or behaviour. 
The basic modelling elements of EF-M are presented in Figure 1. They create a hierarchical product 
architecture representation, which enables the introduction of novel solutions and sub-functions on 
arbitrary levels of granularity into an existing product structure [33]. 

Since EF-M follows the axiom of independence of Suh’s axiomatic design [34], each FR is—in one 
concept—solved by only one DS. If there are more DS per FR modelled, these represent alternatives. 
This is illustrated in Figure 1 in the form of “Solution A” and “Solution B”. 


FR C FR = Functional requirement 
Main function Constraint DS Design solution 
isb C 
icb 


Constraint 





isb_ Is solved by 

rf Requires function 
icb Is constrained by 
iw Interacts with 


FR 
Function #2 
isb 


Figure 1. Enhanced Function-Means (EF-M) modelling elements after [32]. 
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Since new DS are captured for different functions on individual branches of the EF-M tree, they can 
be combined into a multitude of different concepts. In the presented example in Figure 1, only two 
concepts can be instantiated, one with each alternative DS. In Figure 2 two times two alternatives 
become four concepts. 


Product Project EF-M tree Concept 1: Concept 2: Concept 3: Concept 4: 
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Rods in vanes Profile 2 
Figure 2. Instantiation in EF-M and the associated computer aided design (CAD) model. Left hand side 
shows the EF-M model with alternative design solutions (DS), and below the respectively associated 
user defined feature (UDF). To the right are the four resulting concept instances, on top in EF-M and 
their respective embodiment below. All models have been created using a proof-of-concept tool of 
Function and Geometry Exploration (FGE) [16]. 


2.1. The FGE Approach 


The function and geometry exploration (FGE) approach combines the representation of product 
architecture and DR, and the flexibility to introduce novel solutions into an existing legacy design 
with a feature based CAD approach. The approach builds on the function modelling method EF-M 
and combines it with a design automation (DA) approach based on user defined features (UDFs). 
One main assumption behind the approach is that every geometrical element in a product has a specific 
function [35]. This function can be identified, and the related geometry can be isolated. As a result, 
these two information sets can be coupled. The other main assumption for the approach is that if the 
function of a geometry element is identified, alternative solutions for that function can be identified 
and integrated into the product architecture [25,36]. 

The FGE method is based on three main phases: decomposition of the legacy CAD model, 
an innovation stage in the functional domain and an embodiment phase to generate CAD models 
of the concepts. Figure 3 presents the steps inside these three main stages: the decomposition takes 
part in both the function- as well as geometry-model, whereas the geometrical modularisation is 
based on the structure of the EF-M model and the resulting UDF are directly linked to their respective 
DS. The innovation stage includes redesign, where novel functions and solutions are captured in the 
function model, embodiment of said solutions into UDF and the matching of their interfaces towards 
the existing product structure. Lastly, the different concepts are instantiated first from the function 
model and then embodied through the included assembly algorithm. 

To represent the relations between function modelling objects and design features, the FGE 
approach relies on the object model for function and geometry (OMFG). This object model has been 
presented and verified in [16]. 
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Figure 3. IDEFO process model of the applied FGE approach. 


first step the model is decomposed into an EF-M model, in a bottom-up approach where the functions 
of individual geometric elements are identified. Once the geometrical elements are captured as DS, 
they are linked to their function in the form of FR. Potential constraints are also captured in the EF-M 
model in relation to the respective DS. The decomposition of a legacy design into an EF-M model and 
the constraining of the design space in both the functional and geometrical domain is described in 
detail in [37]. Once they are associated to a DS, the geometrical features are captured in a UDF with 
clearly defined interfaces and design parameters. Both of these are captured in the OMFG and provide 
the basis for the DA of FGE. 

Since UDF are available in most commercial CAD systems, these geometry modules can be 
generated and edited by CAD engineers with minimal additional training. Through the flexibility of 
UDF and the possibility to parameterise them, this enables a fine granular control of the geometry via 
the EF-M model [38]. The interfaces between the UDF—and therefor DS—are captured in the OMFG 
and represented in the EF-M model as interacts with (iw) connections. 

Alternative concepts are instantiated from the alternative solutions captured in the EF-M model. 
The instantiation algorithm produces concepts as a combinatorial of all alternative DS, as is illustrated 
in Figure 2. The assembly algorithm based on architectural and geometrical data captured in the 
OMFG automatically generates the CAD models of all technically feasible concepts captured in the 
EF-M model. 

This assembly algorithm, presented in as a UML diagram in Appendix A, Figure A1, 
iterates through all DS of a concept and tries to place the UDF into a part file. While the current 
implementation of the OMFG into a DA tool generates part files, the generation of assemblies is 
theoretically feasible using the same technologies and approaches. The individual UDF are attempted 
to be placed in an arbitrary order, and placed at the back of the queue if the requirements for placement 
(availability of interfaces and parameters) are not given. Only if all UDF have been placed, or all 
remaining UDF have failed to be placed, is the model returned. While this approach is computationally 
expensive, since it may require multiple placement attempts for even larger UDF, it provides a high 
flexibility and theoretical robustness against interface mismatches. 

Following this approach, presented and verified in [16], the geometry not for each new concept 
has to be remodelled, but only for each new DS, which reduces the cost of introducing new product 
concepts based on a legacy design. This modularisation, including a precise interface capture and 
management through the FGE object model, allows for automated generation of the CAD models of 
all feasible concepts. 


3. Materials and Methods 


This study is conducted in tight collaboration with a Swedish aerospace manufacturer, and is 
centred around the design of a new part for turbofan engines, a guide vane (GV). The part is located 
in the bypass of the turbofan engine, positioned just behind the fan, as shown in Figure 4. The set of all 
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GV has the main function to deswirl the airflow from the main rotor, thereby reducing aerodynamic 
losses. Furthermore, the vane has a structural function in that it connects the shroud of the bypass to 
the engine core, thereby creating a load path for the turbine mass to the pylon. Furthermore, the part 
has to withstand axial loads from the engine thrust. 

For the GV being a relatively large component, it is necessary to find low weight designs. As a 
result, low-weigth designs are given high priorities by the design team. This can be achieved by the 
introduction of new manufacturing and material options. The common manufacturing and material 
choice is a metallic core structure [39]. Alternative options are composites or combined solutions. 





l 
<5 


Figure 4. A Rolls Royce Trent 900 engine rendered in CAD with one GV highlighted in orange 
(CAD model by Chris Shakal on grabcad.com). 


A GV commonly consists of a vane structure, a predominantly aerodynamically defined structure, 
and two end attachments which fasten the vane to the engine core on the inside, and the outer fan case 
on the outside, as can be seen in Figures 4 and 5. Both parts provide interfaces the GV assembly has to 
comply to. 

The development team is composed of both design and analysis engineers with experience in 
aerospace engineering ranging from five to thirty years. These engineers partake in the development 
in defined percentage of their work-time, between 10% and 70%. 

It is important to mention that the company is using and developing KBE techniques with a high 
degree of DA for design space exploration. This means that the engineers possess a benchmark system 
for comparison with the FGE. The benchmark system can be compared in make-up and extend to the 
parameterisation-based multi-disciplinary analysis (MDA) framework presented in [40]. 


Study Setup 


The research has been performed as a participatory action research study, where the researchers 
accompany the development process investigating into new opportunities for a fan frame GV [41]. 
The researchers participated in regular team meetings together with the design team, and had the 
ability to observe the engineers during their work. Notes from these meetings and observations are 
part of the data supporting the presented results. 

Beyond the observation throughout the product development project, a set of studies has been 
conducted to investigate the perception of the presented method. For the specific purpose of validating 
the FGE method and tool, two consecutive workshops were organised together with the design team. 

In a first workshop together with the GV development team (eight participants), a functional 
decomposition of the part was performed and captured. Hence, this workshop is referred to as the 
“decomposition workshop”. The actual workshop was predated with a one hour long test-workshop in 
one of the regular team meetings to evaluate the motivation of the practitioners and the requirements 
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for the actual workshop. The results of this “test-workshop” are included in the results of the 
decomposition workshop. The decomposition was performed using large-scale print-outs of drawings 
a GV which were annotated in a group exercise. The annotated drawings were collected, digitised and 
used for the creation of function-geometry model in FGE for the legacy design. 
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Figure 5. Example of function and design solutions of a guide vane (GV), connected to geometry 





elements as UDF. The different geometry elements are colour coded to show the extend of the UDF. 


A second workshop (eight participants) covered the innovation stage of the FGE approach. 
The workshop was held with the same development team, using the same product. The workshop 
was held online due to the social distancing regulations in response to the COVID-19 pandemic. 

For the second workshop, a proof-of-concept tool of the FGE approach was created. The tool 
works server-based with a web interface (illustrated in Figure 6), allowing users with different skill 
levels and operating systems to connect to the same database. The tool incorporates a fully functional 
implementation of an EF-M modeller with additional functionality to capture geometric information 
for each DS. As such, it is able to capture alternative design solutions and instantiate them into 
different concepts. Furthermore, the tool includes a geometrical modelling algorithm to automatically 
generate CAD models of each concept, based on previous geometrical modularisation in Siemens NX. 
The structure and mechanics of the tool are illustrated in [16]. 

After each workshop, the practitioners were asked for their experiences through questionnaires 
and open feedback discussions. The results of the first workshop were recorded through adhesive 
notes and remarks on the drawings, as well as protocol notes. The results of the online workshop were 
captured through audio and video recording and through a change protocol in the database of the 
FGE tool. 

For this publication, a generic design of a GV was created, which differs from the design used 
in the workshops. However, the generic design has a sufficient degree of realism to enable dialogue, 
but with company specific features left out. 
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Figure 6. Web-based interface of the FGE tool “omfgDSE” (object model for function and geometry 
based design space exploration), showing a simplified EF-M model of a GV. The right-hand pane 
illustrates the design rationale of the selected DS, such as associated function, geometry and parameters. 


4. Results 


This section presents the results of the application of FGE in the industrial context as well as the 
practitioners feedback and observations. 


4.1. Design Space Exploration Using FGE 


The DSE approach is separated into three phases decomposition, innovation and embodiment, 
following [36]. The presented workshop focuses on the interactions with practitioners in the phases 
decomposition and innovation. 

In the first workshop (decomposition), the legacy design of the GV was captured in an EF-M 
model based on the available geometry. The model contains three main functions, and 18 different 
design solutions on the concrete (lowest) level. Since the model was created through decomposition, 
all DSs are directly linked to one or several geometry elements. A modularised CAD model was 
created based on the original geometry, modularised according to the decomposition. 

Figure 5 exemplary shows this modularisation and connection for four individual DSs and the 
related UDFs. For example, the leading edge of the decomposed GV is covered by a metallic edge. 
This metallic edge, as a geometrical element, has been associated with the function “Withstand impacts” 
through the solution “Metallic edge”. The respective geometry, of both the edge insert and the cutout 
in the vane to place it in, are captured as two UDFs which are then associated with the DSs. 

This decomposition served as input for the following workshop, focused on innovation. 
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In the second workshop (innovation) the prototype tool was used to capture novel design solutions. 
As a result, the product model was extended by five new sub-functions with their respective solutions 
and 10 new alternative DSs. As an example for new FR, sub-functions related to acoustic performance 
were added to the model. The nature of the DSs ranged from material alternatives over different 
interface approaches to variations in shape and dimensions of the existing solutions. The entire EF-M 
tree, with all new FR and DS highlighted with red borders, is shown in Figure 7a. Based on the new 
alternative DSs, 52 different concepts could be instantiated. 

The majority of the novel solutions shows a conceptual difference from the legacy design. As an 
example, the FR “Join vane to attachment” is in the legacy design solved by the DS “Adhesive 
connection”, meaning that the vane is glued to the attachment using resin. This DS, and the respective 
UDE are illustrated in Figure 5. The two alternative DSs, as results from the workshop, were “Bolted 
joint” and “Fully integrated solution”. The respective DS are shown in Figure 7b, whereas the geometry 
related to these DSs is shown in the context of the GV in Figure 8. 

The embodiment phase of the FGE approach has not been performed in its entirety, since the DA 
section of the proof-of-concept tool has not sufficiently matured to handle the complex geometries of 
the solutions devised in the innovation workshop. For a verification of the proof-of-concept tool on 
another turbine structure, see [16]. 

Figure 8 shows an exemplary embodiment of three concepts with alternative DSs for the FR 
“Join vane to attachment” shown in Figure 7b. While otherwise using the same configuration, 
each of the geometries is adapted to accommodate the respective geometric changes in the concept. 
The instantiation and assembly algorithms embedded in the FGE model allow for automated generation 
of the CAD models of all feasible concepts. 








Configurable component 











(a) (b) 
Figure 7. EF-M model after the innovation workshop. (a) Full EF-M tree, with all new DS and functional 
requirement (FR) highlighted with a red border. Readability has been dismissed on purpose to protect 
company intellectual property (IP). (b) Detail: alternative solutions for DS “Join vane to attachment”, 
highlighted in blue in (a). 














Figure 8. Embodiment of three concepts, differentiating in the solutions for the FR “Join vane to 
attachment”, from left: (a) “Adhesive connection”, (b) “Bolted joint” and (c) “Fully integrated solution”. 
Respective alternative geometric elements relating to the respective user defined feature (UDF) are 
highlighted in colour. 


The following section reports on the practitioners’ feedback and reflections after the application 
of the method on their industrial use case. 
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4.2. Reflections Form Industrial Practitioners Regarding the Application of the Method 


The practitioners have stated that the workshops—as activities in addition to their regular 
design activities—contribute to the development project. Furthermore, all workshop participants 
would like to repeat the workshop concept, both in the same, but even more so in other projects. 
The following sections will go through the main mechanisms that have been highlighted by the 
industrial practitioners, in particular the ability of the approach to: 


1. Represent and capture design rationale and the design space 
Capture, integrate and model novel solutions 
3. Provide support for the embodiment of novel concepts that would otherwise remain unexplored 


4.2.1. Capture and Representation of Design Rationale and Design Space 
“T learned a lot of things about the GV just by reading the graph.” 


—aerodynamics analysis engineer 


This quote illustrates a common problem during the design of aerospace products. As components 
are designed and conceptualised over long times, and by heterogenuous teams, it is difficult to retain 
the design rationale (DR) of a product (that is, why a product has been designed in this form) [42]. 
This DR is especially hard to read in highly integrated and complex products such as aerospace 
engine components, where a single geometrical artefact may fulfil multiple functions in different 
domains [43]. For practitioners, the FGE approach supports the ability to visualise the rationale in both 
the geometrical domain, visualising “what is”, and the functional domain, “why it is”. 

In terms of the use of FGE, it appears to be a major benefit of FGE for the practitioners that 
they improved their knowledge about the product and participated in a forum that permitted 
exchange about the product’s functions and solutions. The questionnaire accompanying the workshop 
revealed that no team member had full confidence about the GV, but most rated their understanding 
of the GV higher after the workshop than before. The workshop also fostered the awareness of 
meta-product-knowledge, since many engineers were not even aware who knew the rationale behind 
a specific element. Only through discussions and participation of engineers of multiple disciplines 
could the entire DR be captured in the EF-M model. 

An example for how FGE enables access to DR is the interface of the tool used in the study, 
which collects all information of a solution in a single place. Figure 6 shows on the right hand side how 
the geometry, parameterisation and other information about a solution is collected in a single window 
pane. The top of the pane shows which function is solved, and the highlighting of the solution shows 
its position in the product architecture. 

Furthermore, the approach can illustrate the geometrical interfaces between solutions through 
iw connectors in the EF-M model. This are seen as black arcs in the interface of Figure 6, 
left hand side. This approach of information representation was referred to as “comprehensive” 
by the workshop participants. 

Beyond the information about individual solutions, the FGE interface provides an overview over 
the possible alternatives through EF-M. This was perceived to be an efficient way to represent the 
design space, providing “good possibilities to look at different solutions”. However, it was also noted 
by one participant in the workshop that the function models can quickly grow to be too complex. 
The function model created in this study (see Figure 7a) was considered “almost too large”. 
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4.2.2. Ability to Capture, Integrate and Model Novel Solutions 


“Parametric solutions exist but limited to parameters that we easily can change, dimensions, 
angles or thicknesses for example. When it comes to exploring different design solutions in 
early phases we foresee larger challenges.” 


—development team leader 


This quotes stands exemplary for the practitioners’ reflections over the in-house KBE system for 
DSE adopted by the company, which is built around state of the art KBE approaches. However, 
the practitioners are well aware about the limitations of highly automated CAD approaches in 
representing conceptually different design solutions. While the exploration of dimensional variety 
is possible and appreciated by the practitioners, the need to realise topologically and geometrically 
different options has been stated explicitely. This matches with observations from literature, where it 
is stated that “a parametric model offers little flexibility” [44]. 

According to the practitioners, FGE has the potential to deliver this through the modularisation 
of the geometry, and the connection to the DR of the product. 

It was perceived by the participants of the workshop that the use of the FGE approach supports 
the introduction of new solutions (4 on a scale of 1 “no support for new solutions” to 5 “totally 
supports”). It was furthermore stated that part of the identified DS would probably not have been 
found—or considered—in a regular modelling approach, with a majority (71%) stating “maybe some, 
but not all”. The other participants stated either “none of them” or “most of them” would have been 
identified without FGE. 

Beyond only exploring new solutions for existing functions, FGE also enables the capture of 
newly identified functions or sub-functions. The workshop results showed that through five new 
functions being identified during the innovation workshop. Practitioners feedback further supported 
this, although the ability to integrate new solutions is perceived to be more important than that of 
new functions. 

To summarise this, the practitioners perceive FGE as a viable method to capture novel solutions 
and functions. Having a formal, and even computer supported, modelling approach at hand when 
generating conceptual alternatives directly in creative design meetings was found very applicable. 


4.2.3. Support for the Embodiment of Novel Concepts 


“There are still gaps within early concept phase—this is one possibility to generate, 
and evaluate, lots of concepts, get trend curves and so on.” 


—design engineer 


As has been stated in the section above, there is a substantial need for a DSE modelling approach to 
capture geometrically and topologically different DS. The practitioners agree that FGE is an approach 
which can support such exploration. 

However, the sole representation of alternative concepts in function models would be insufficient. 
To enable the necessary analyses, geometry models in the form of CAD are required. 

The practitioners stated that the CAD models needed for effective DSE need to be of sufficient 
quality to be meshed automatically. This requires certain quality criteria such as surface continuity. 
Furthermore would the CAD models require automated “tagging” of relevant geometry elements for 
analyses, for example where forces attack or constraints apply. 

After a demonstration of the assembly algorithm of FGE and the geometry modularisation, 
the practitioners agreed that such an approach can deliver the required models. Therefore, the ability 
to generate topologically different variant concepts was seen as a benefit of the presented method: 


“Generation of CAD models based on different configurations is a key functionality.” 


—development team leader 
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4.2.4. Implementation Challenges 


Despite the interest in the approach, the practitioners pointed out certain implementation challenges. 

As a feedback to the modelling process, using EF-M to capture and represent alternative solutions, 
a practitioner mentioned that the re-use of existing solution’s geometry for another function is 
hardly possible to model using EF-M. Although there exist horizontal relations in EF-M modelling, 
they increase the modelling complexity and reduce readability. Therefore, such complex relationship 
modelling has not been implemented yet in the prototype tool, and does not contribute to alleviate the 
challenges of complexity in function modelling stated by [31,33]. 

The ratio of complexity between the functional domain and geometrical domain was point for 
critique. It was stated that a large function model (such as the one in the workshop became through 
the innovation) reduces the overview over the model. It was suggested to use less FR/DR elements 
and more, but potentially more complex, geometric modules. However, the experience with function 
modelling and/or web-based applications may play a role in how difficult the tool use was perceived. 

While FGE supports the exploration of technically, functionally and topologically different 
solutions, the in-house KBE solution is optimised for parametric variation. Being used to the already 
existing MDA method at the company, the practitioners rather viewed an implementation of FGE as a 
complement, preferably integrated with the existing system. 


5. Discussion 


The FGE approach aims to support DSE in the conceptual product development phases through 
a three-phase approach, namely decomposition, innovation and embodiment. The approach relies 
on the OMFG, which provides a connection between the previously divided modelling domains of 
function and geometry. As a result, novel solutions can be integrated into an existing legacy design, 
and be embodied based on information stored in the function model. This information, the DR of the 
product, is the driving force behind the implemented DA approach. The coupling of the function and 
geometric modelling domain, and the resulting opportunities for DSE, distinguish the approach from 
other DA approaches and CAD model parameterisation. The approach has been evaluated in a case 
study with a product development team of an aerospace manufacturer. 

The practitioners gave an overall positive feedback about the approach, stating that FGE—or 
a similar approach—can provide what is needed to improve current product development practice. 
The need was stated to be about capturing and presenting alternative solutions in a systematic way, 
as well as an automated generation of CAD models ready for analysis. The use of a non-geometrical 
modelling approach representing teleological product knowledge, “how it all connects”, has been 
received positively by the practitioners. This is in line with authors such as [8,14,15] who call for 
models for DSE to carry information beyond the geometrical domain. 

While the results of this study can be interpreted as generally positive for the presented FGE 
method in terms of user acceptance and improved DSE, it is only a single case study, and hence difficult 
to generalise from [45]. While the engineering design research community is aware of the challenge 
to properly validate a new method or methodology [46,47], this study can only be an initial step into 
validating FGE. Similar studies with other development teams, and in other engineering domains, 
are required before a conclusive image of the contribution of FGE to DSE can be presented. 

Practitioner feedback has shown challenges of understanding the instantation, EF-M and CAD 
integration of the approach. In general, from the observations throughout the workshops, did the 
time to understand and learn the new approach take up a large part of the workshops’ time. Hence, 
the studies results and conclusions are influenced by the novelty of the tool and method for the 
practitioners. This can be a factor in the enthusiasm and positive feedback for the method, in a variant 
of the “Hawthorne-effect” [48]. However, it could also be argued that the method can only show 
its full contribution once the practitioners have grown familiar with it. Which of these two factors 
predominates cannot be said at this moment. 
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Many participants highlighted the contribution of the method for knowledge capture and 
exchange inside the group. Although this is not among the primary objectives of the FGE method, 
the capture, representation and retrieval of design rationale is beneficial to the product development 
process [49]. However, it can be argued that not the use of this specific method, but the general 
interaction of engineers of different disciplines focused on a shared product could provide this [50]. 

The registered numbers for captured novel design solutions, concepts and functions have 
been identified by the participants in the study as subjectively higher than without the method. 
However, to be able to judge a significant impact of the use of method on the quantity—or 
quality—of the generated design concepts, further and especially more longitudinal studies are 
required [45]. The presented FGE approach for DSE relies heavily on a DA implementation for 
the assembly of individual concepts’ CAD models. This approach has been presented in detail 
in [16]. The practitioners in the study have already voiced concerns about the robustness of the DA 
approach and the quality of the resulting CAD models for meshing and FEM. Similar to other DA 
approaches and parameterisation, FGE is sensitive to modelling errors, ambiguously defined datums 
or over-constrained models [20,21,23]. The definite sensitivity to modelling errors has to be evaluated 
in further studies. 

However, the actual contribution of the FGE approach lies not in the DA implementation, which in 
similar fashion can be found in existing feature-based design approaches [51]. More than providing 
a feature-based geometry approach does FGE provide a function-based approach to geometry 
modelling. A change in product function—or change of a solution—directly impacts the CAD model. 
The underlying OMFG provides a direct link between the two domains, and therefor enables a 
DSE approach which is guided by “the most important concept in determining a product’s basic 
characteristics”: function [14]. 


6. Conclusions 


FGE is a novel DSE approach which has been presented in [16]. While this previous presentation 
has presented the technical workings of the approach and verified the DA algorithm, the approach as 
a method in EDR has to be validated through application in an industrial context. 

This paper presented the application of a DSE approach combining function modelling and 
automated geometry model generation for the conceptual product development phases of turbofan 
engines. The modelling approach tackles the challenges of representing alternative concepts in the 
design space, design rationale capture and storage as well as support for the generation of CAD 
models. Through the application of the FGE approach in an engineering case, the method has been 
validated to target a real problem in the conceptual development of aerospace engine components and 
to provide a feasible solution for it. 

Through the use of a prototype tool implementing the FGE method, engineers in a aerospace 
supplier were able to share product knowledge and explore novel design solutions. This has lead to 
both, a perceived increased knowledge about the developed product inside the development team, 
as well as the opportunity to explore previously unconsidered design solutions. According to the 
statements of the product developers participating in this study, this leads to more and more novel 
functions and solutions being explored. 

The FGE approach uses function modelling as a tool to represent the design space and perform its 
population. While this has shown to be feasible, domain specific challenges similar to as reported for 
example in [31] have shown. This points towards further development for EF-M modelling in order to 
increase applicability and reduce abstraction. Furthermore did practitioners point out the problems of 
modelling a “make use of” relationship for a DS, as described in Section 4.2. Similar challenges have 
already been identified in [33] and pointed towards further development of the EF-M approach. 

The presented prototype tool has shown to efficiently provide an interface for collaborative DSE. 
To a degree it has reduced the abstraction of function modelling by providing a connection to tangible 
geometry. However, further studies with increased practitioner training are seen as a goal for future 
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work. While the embodiment process has been illustrated in [16], further studies have to show the 
scalability of it. In particular, the quality of the generated CAD models—with a focus on meshing—has 
to be in the focus of development, since in the presented study it has been highlighted that it is crucial 
for analysis of the concepts. 

In summary, while there have been identified areas for improvement, the FGE approach was 
perceived to be “one possibility to generate, and evaluate, lots of concepts”—which is what it aims for. 
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Abbreviations 


The following abbreviations are used in this manuscript: 


CAD computer aided design 

DA design automation 

DR design rationale 

DS design solution 

DSE design space exploration 

EDR engineering design research 

EF-M enhanced function-means 

FGE function and geometry exploration 
FR functional requirement 

GV guide vane 

iw interacts with 

KBE knowledge-based engineering 
MDA — multi-disciplinary analysis 

OMFG _ object model for function and geometry 
UDF user defined feature 
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Appendix A. UML Process Model of FGE Assembly Algorithm 
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Abstract: A comparative sensitivity study for the flutter instability of aircraft wings in subsonic 
flow is presented, using analytical models and numerical tools with different multidisciplinary 
approaches. The analyses build on previous elegant works and encompass parametric variations of 
aero-structural properties, quantifying their effect on the aeroelastic stability boundary. Differences in 
the multifidelity results are critically assessed from both theoretical and computational perspectives, 
in view of possible practical applications within airplane preliminary design and optimisation. 
A robust hybrid strategy is then recommended, wherein the flutter boundary is obtained using 
a higher-fidelity approach while the flutter sensitivity is computed adopting a lower-fidelity approach. 


Keywords: aircraft design; aeroelastic stability; sensitivity analysis; flexible wing; subsonic flow 


1. Introduction 


Within aircraft multidisciplinary design and optimisation (MDO) [1,2], efficient methods and 
robust tools are highly sought after as effective reduced-order models (ROMs) [3-5] for fast parametric 
aeroelastic analyses [6-8], especially for sensitivity and uncertainty evaluation purposes at the 
preliminary design stage [9-11]. Simplified semi-analytical formulations for the bending-torsion 
instabilities of flexible wings have been available for a long time but have inherent limitations [12-14], 
whereas complex fluid-structure interaction (FSI) [15] procedures based on the finite element method 
(FEM) [16] and computational fluid dynamics (CFD) [17] have recently been used to improve accuracy 
and generality, but remain computationally expensive and require special care [18,19] with significant 
effort to pre- and post-process the numerical simulations. 

Continuing previous studies on multilevel techniques for practical aeronautical applications [20,21], 
this conceptual work investigates a robust hybrid strategy where the flutter boundary is accurately 
obtained using a higher-fidelity approach while the flutter sensitivity is efficiently computed by 
adopting a (tuned) lower-fidelity approach as an effective ROM. A comparative sensitivity study is 
hence presented for the aeroelastic stability boundary of a uniform cantilever wing in subsonic flow 
as the standard benchmark [22]. The Typical Section idealisation [23,24] is employed as the analytical 
lower-fidelity model, whereas the numerical higher-fidelity model couples a beam-based FEM with 
panel-based CFD. Considering Loring’s wing [25] with a NACA0002 airfoil [26], the flutter analyses 
build on earlier publications [27,28] and encompass parametric variations of wing properties, such as 
its structural inertia, stiffness and aspect ratio. The effects of the latter on the aeroelastic instability are 
quantified and differences between lower-fidelity and higher-fidelity results are critically assessed from 
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both theoretical and computational perspectives, in order to study the necessary trade-off between 
complexity and costs of model fidelity for intensive industrial MDO activities [29]. 


2. Problem Formulation 


Uniform cantilever wings have long been used as the standard benchmark for fundamental 
studies on aeroelastic analyses and parametric sensitivities for the bending—torsion instabilities of 
flexible wings, both numerically and experimentally [12,22]. They feature constant material properties 
(i.e., structural density 0m, Young’s elastic modulus E and shear modulus G), chord c, mass m and 
moment of inertia p per unit length at the inertial axis xcg (i.e., the line of sectional centres of gravity); 
bending stiffness EI and torsion stiffness GJ at the elastic axis xp, = 0 (ie., the line of sectional shear 
centers) along the semi-span /; structural damping is typically small and conservatively neglected [30]. 

With 6(y,t) and h(y,t) the pitch and plunge motion of the flexural axis at the time t, 
respectively, the wing vertical displacement w(x, y,t) is given as w = h+ x0 directly, x and y being 
the chordwise and spanwise directions. Assuming an Euler—Bernoulli beam model [31], the linear 
system of coupled partial differential equations (PDEs) for the dynamic aeroelastic response of wing 
bending and torsion undergoing small deformations reads: 


se i se 


m (h _ xcc6) +FEIh = AL, ud — MXCG (h _ xcc6) — GJh = AM, (1) 


with AL(y,t) and AM(y,t) being the sectional lift and pitching moment, respectively; the governing 
equations are then completed by both geometrical and natural boundary conditions, namely: 


6 (0,t) = GJ@ (1,t) =0, h(0,t) =h (0,t) =ElIh (1,t) =EIh (1,t) =0. (2) 


As for the case of Loring’s wing, the problem formulation assumes an isotropic material without 
loss of generality; an anisotropic material (e.g., for the case of composite wings [32]) would feature 
elastic coupling between bending and torsion (due to different mechanical characteristics in respective 
directions [12]) but would not introduce conceptual changes to the overall methodology and main 
conclusion of this work. 


2.1. Uncoupled Natural Vibration Modes 


Assuming h © ¢,x;, and 0 © dgx9, where x(w,t) are the generalised coordinates of the uncoupled 
vibration mode shapes ¢(y) with natural frequencies w, the latter for the cantilever wing’s bending 
and torsion dynamics can independently be obtained by solving the relative homogeneous PDEs via 
separation of variables with the respective boundary conditions, namely [12]: 


2 (Ed 
(o4.= (+) —, On = . el coshycosy+1=0, cosv = 0, (3) 


m 
py = cosh (1¥) ~cos(o4) — (SERT*S™T) [san (y¥) —sin(o¥)],— go=sin(o¥), 


which form a complete set of modal bases for the generalised solution of the aeroelastic equations; 
yet, note that these uncoupled bending and torsion modes are inherently orthogonal within their own 
bases but not to one another. 


2.2. Unsteady Aerodynamic Model 


According to thin aerofoil theory for incompressible potential flow, the sectional unsteady 
aerodynamic loads due to the wing motion may be written as [22]: 


A= =P Ez (U6 — hh + xycb) + KUC)/_ (Crv) » “VO peg. (5) 
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2 
AM = —56 154 Be + xcpUé — XMC (h _ xmc) a KUX AcC] /q (crv) | , (6) 


where p¢ and U are the reference air density and horizontal flow speed; the aerodynamic centre x 4c 
(where the circulatory load acts) and the control point xcp (where the non-penetration boundary 
condition for the vertical flow speed V(y,t) is imposed) lay at the first- and third-quarter chords, 
respectively; xc refers to the mid-chord (where the non-circulatory load acts); Cj;, is the derivative 
of the sectional lift coefficient C; with respect to the nominal angle of attack a. 

Within tuned strip theory [21], the scaling factor « introduces steady three-dimensional downwash 
effects due to the tip vortices [33] and takes part in the derivative of the wing lift coefficient 
Cr /y = KC/, [34], while the complex lift-deficiency function Cr(k) is defined in the reduced-frequency 
k domain and embeds the lag due to two-dimensional inflow effects from the travelling flat 
wake [35], namely: 


4e 7TH 
C — 270 1 — " ) J k= y C 
ae ( 3/3 7tne + Ch/q : 


with e€ being the aerofoil thickness ratio, 7 the wing aspect ratio, e(7) Jones’ edge-velocity factor (i.e., 


H? cw 
= —- _ k = ey 7 
Hy ing 2U ”) 


the ratio between wing semi-perimeter and span [36]) and H?(k) Hankel’s functions of the second 
type and n-th order. Note that the proposed approximate correction for the wing lift coefficient is 
analogous to applying Oswald’s efficiency factor [37], whereas x = 1 for standard strip theory. 

In summary, unsteady aerodynamics account for both circulatory lag and non-circulatory 
inertia effects [38], whereas quasi-unsteady aerodynamics neglect circulatory lag terms [39] within 
a quasi-static approach (i.e., Cr ~ 1 with k << 1). Quasi-steady and steady aerodynamics then 
disregard non-circulatory inertia terms too, within a static approach (i.e., Crp = 1 and k = 0) where 
the control point coincides with the elastic axis in the former case (in order to avoid unrealistic 
aerodynamic damping [23]) and all time derivatives are eventually discarded in the latter case. Note 
that the structural motion represents a time-dependent boundary condition for the incompressible 
airflow in all cases, but a synchronous proportional variation of the circulatory airload occurs within 
quasi-static approaches only [40]. 


2.3. Modal Approach and Stability Analysis 


Given the natural vibration modes of the wing, the relative equations of motion can be transformed 
into ordinary differential equations (ODEs) by employing Ritz’s method [41], where generalised 
coordinates {x(t)} multiply the natural vibration mode shapes ¢. Recasting the whole aero-structural 
model in vector-matrix form, the governing modal equations then read: 


IM} ix + IC) ixs + [K] tas = 10, {Fa} = [Mal ixt + [Calixt + [Kal txt, (8) 
where the aeroelastic mass |M], damping |C] and stiffness [K] matrices are given by: 
[M] = [Ms] — [Mal], [C] = [Cs] — [Cal], [K] = [Ks] — [Ka], (9) 


and enforce a monolithic FSI, with {F,(t)} the aerodynamic load; as anticipated, [Cs] = [0] is hereby 
assumed for convenience without loss of generality. 

Due to linearity, the parametric aeroelastic stability of the wing is then monitored from the root 
locus of the characteristic equation for the flutter determinant as the reference airspeed increases; 
in particular, the system becomes metastable when the real part of at least one eigenvalue A vanishes 
and leaves the response undamped, namely: 


fx} {up eM, (2? [M] + A [C] + [K]) {u} = {0}, (10) 
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with A = iw in the case of dynamic flutter where a coupled resonant harmonic motion excites the 
eigenmode {u} (Hopf’s bifurcation occurring when a couple of complex conjugate eigenvalues cross 
the imaginary axis), whereas A = 0 in the case of static divergence (occurring when the aeroelastic 
stiffness matrix becomes singular). Note that the eigenvalues and eigenvectors are suitably assumed 
to form a complete distinct set, as per most of the practical cases [42]; of course, a sufficient number 
of natural vibration modes shall be employed in order to grant proper convergence of the modal 
analysis [43]. 


2.4. Sensitivity Analysis 


The aeroelastic stability boundary acts as a typical constraint in aircraft MDO problems, where 
parametric sensitivities of both flutter and divergence speeds with respect to variations in the design 
variables are required for gradient-based optimisation algorithms [44,45]. In general, such derivatives 
may numerically be obtained by performing many simulations with small alterations of each and every 
design variable; nevertheless, the aircraft’s conceptual and preliminary design stages would require 
extensive robust computations for the trends of both objective function and constraints while changing 
the design parameters [46,47]. Effective analytical solutions (with explicit expressions in a few special 
cases) have hence been derived by differentiating the governing equations, taking advantage of the 
self-adjoint nature of complex eigenproblems and then separating real and imaginary parts [48,49]. 

Especially when the sensitivity of the critical mode shape is also required, the continuation method 
may efficiently be employed to calculate the critical point and obtain all desired derivatives at the 
same time [50]. The eigenproblem arising in flutter and divergence analysis may indeed be considered 
as a system of nonlinear equations with a normalisation condition for the eigenvectors; the solution 
is then numerically sought by the Newton—Raphson method, where the aeroelastic equations are 
differentiated with respect to a design parameter of interest p and the system is linearised at each 
iteration [51]. Within a single unified procedure, the eigenvalues’ and (right) eigenvectors’ derivatives 
may then automatically be determined along with the eigensolution itself at once, with no information 
about the transposed (left) problem being required; higher-order derivatives may be obtained by 
further differentiation. In particular, with U = U,(p),w = w-(p) and {u} = {u-(p)} in the critical 
flight condition, the nonlinear eigenproblem for the flutter mode reads [9]: 


[S]{ue}=0, {uc} {ue} = 1, [S] = —w; [M] + iw [C] + [K], (11) 


and is differentiated to give the eigensolution sensitivity with respect to the design parameter as: 


Re [S| —Im|S] —Re{v} Re{c} Re {uc} Re {1} 

Im [S| Re [S| Im{v} Im{c}} 0 ) Im{uc} _ J Im{ey (12) 
Ref{uc}" —Im{uc}" 0 0 op U, 0 ‘ 
Im{uc}' Re {uc}* 0 0 We 0 


where the submatrices and subvectors involve the derivatives of the aeroelastic model, namely: 











icr= mn {ets al = -w2 + iW, < + “ — 2w, |M] +i[C], (13) 
(y= ud, FE = wD i D+ (ug) 
i; = -> {uc}, = = -w3 Fee + oI (15) 


in the degenerate case of static divergence with Im {u-} = 0 and w, = 0, all quantities are inherently 
real and the matrix of algebraic equations above specialises without the second and fourth rows and 
columns, respectively. Following previous studies [27], the derivatives may finally be normalised with 
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respect to the reference values of the involved quantities, thereby comparing all sensitivities within 
a coherent representation. 


3. Lower-Fidelity Model 


The Typical Section abstraction [52] is hereby employed as the analytical lower-fidelity model, 
as it is conceptually illuminating, inherently robust and computationally efficient within its limitations. 
It resembles the section of the uniform wing and its fundamental aeroelastic behaviour for the 
bending-torsion coupling mechanism; matching the wing inertia and natural frequencies, such an 
abstraction, may be regarded as a condensed ROM. The structural arrangement is a classic mass-spring 
system for the wing section, where the coupled pitch and plunge motion is restrained by vertical 
and angular linear springs with equivalent stiffnesses k, = mw? and kg = pws (see Appendix A), 
respectively; it is representative of the inertial and elastic properties per unit length of the wing at about 
75% of its span, where the structural arrangement becomes progressively flexible and the aerodynamic 
load is still high. The model is hereby extended to include all bending modes having natural frequency 
below that of the first torsion mode, as the instability may involve any combination of them in general. 

Considering the first two uncoupled bending modes and the first uncoupled torsion mode of 
Euler—Bernoulli’s beam model (with 7; * 1.875, y2 ~ 4.694 and v, & 1.571) [12], the aeroelastic 
equations of motion for the wing’s Typical Section in (time-varying) steady incompressible flow are: 


m 0 —mxccf kp 0 0 
[Ms] = 0 m —mxccg |, [Ks] = |0 kyr O|, (16) 
—mxcof —mxccg Pw+Mxa 0 O kg 
00 f 
[Ma] = [0], [Ca] = [0], Ka] =5UPcCrja |0 0 g |, 7 
0 O —XAC 


where the cross-projections of the non-orthogonal modes scale the off-diagonal coupling terms: 


1 il 1 /! 5 
= a Gn (¥, 11) 0 (Y,M1) dy, = i n (Y, 2) Po (Y V1) dy, C me (18) 


in particular, all the latter are constant and read f ~ 0.959, g = 0.274 and r & 39.275 for homogeneous 
aero-structural properties. Note that assuming a static aerodynamic approach is consistent with 
the rigorous application of the scaling factor x, which is hereby based on steady lifting-line theory; 
no damping is hence provided to the aeroelastic response. It is also worth stressing that the equations 
for the first and second bending/plunge modes are not coupled directly, but the equation for the 
torsion/ pitch mode couples them indirectly. 


Aero-Structural Parametric Derivatives 


The sensitivity of the aeroelastic matrices with respect to aero-structural parameters can readily 
be obtained in explicit analytical form, where the chain rule applies [28]. In particular, the derivatives 
with respect to the material density 0, and elastic modulus E are given by: 

O([K] 1 


= |0], >F FE [Ks], (19) 





d|M] 1 O[K] | 0 |M] 
Don el bo FE 


whereas those with respect to reference flow speed U and perturbation frequency w read: 


0 [M] 


0 [M] 
Ou Ow 


OW 


=], Se=sKKal, as 


TG = (0), = (01, (20) 
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and with respect to the wing semispan | it is: 








20 0 
a a[K) 2 1 ac 
aM — (ay, AN = —5 [D] [Ks] + —— [Kal Lia [DJ=|0 2 of, =) 
noe 001 
ACL /» _ 2 OC /y ACL /« :: 1—« 
at ey “ay a) mi 


4. Higher-Fidelity Model 


The high-fidelity model is based on a FEM of a slender beam solved with the commercial code 
Nastran [53] for the structural analysis; the latter is coupled with either the doublet lattice method 
(DLM) available in the same software [54] or an in-house panel code based on Morino’s boundary 
element method (BEM) for the steady and unsteady aerodynamic analysis [55,56], thereby enforcing 
a numerical FSI. Both Nastran-based structural FEM and in-house aerodynamic BEM are generated 
with automatised routines, in order to ease parametric studies; this capability is exploited to perform 
numerical convergence studies (see Appendix B) and compare the analytical flutter derivatives with 
their numerical counterparts obtained via finite difference. 


4.1. Structural Model 


Following previous works [21], the wing structure is modelled with CBEAM elements and 
accounts for the distance between inertial and elastic axes; the node lying at the wing root is then 
clamped (see Figure 1). Further FEM nodes are placed at the leading and trailing edges of the wing 
and connected to the beam nodes with the rigid elements RBE2, in order to support mapping with the 
aerodynamic grid (in terms of both structural deformation and aerodynamic load: the latter modifies 
the former, which changes the airflow boundary condition in turn) according to the closely-spaced 
rigid diaphragm (CSRD) assumption [57]. For the natural vibration analysis, shear deformation is 
neglected in order to obtain the Euler—-Bernoulli beam model, with PBEAM defining the properties (i.e., 
inertia and stiffness) of the beam element and SPC1 defining the single-point constraint for the clamped 
root. Using Nastran’s SOL103 to obtain the structural eigenvalues and eigenvectors, the vibration 
analysis is performed while selecting Lanczos’ method [58] as available in EIGRL and normalising the 
modes to unit values of the generalised mass. 





Figure 1. Structural FEM of Loring’s wing. 
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4.2. Aerodynamic Model 


Standard flutter prediction methods in the industrial environment are based on aerodynamic 
panel codes (mostly DLM [59,60]) that idealise wing and empennage as lifting surfaces; if considered, 
the aircraft fuselage is treated with dedicated elements for non-lifting bodies. The lower computational 
cost compensates for the lower fidelity of the model; however, the reliability may be increased by 
correcting the aerodynamic influence coefficients with higher-fidelity data (which are generally not 
available at the preliminary stage of an MDO process though). Considering small perturbations of 
an unsteady irrotational flow, the BEM proposed by Morino [55,56] is hereby adopted and consists of 
an integral representation of the velocity potential at any point of the computational domain in terms of 
the values on the surface surrounding the aircraft body and wake; the principle of superposition applies. 
The linear equation of acoustic waves propagation governs the unsteady flow for the case of isentropic 
compressible fluid and accounts for the finite speed of sound [61], whereas Laplace’s equation is 
resumed for the case of steady flow by taking advantage of Prandtl—Glauert’s transformation [62]. 
Green’s function representing a unit-impulsive point source, the theoretical formulation of the 
boundary value problem for the perturbation potential is then based on Green’s formula [33], 
with (Neumann-type, generally instationary) non-penetration boundary conditions on the aircraft 
surface and wake as well as (Dirichlet-type, generally stationary) unperturbed asymptotic conditions 
far from the latter. Bernoulli’s theorem is linearised to calculate the pressure coefficient and Kutta’s 
condition is imposed at the trailing edge of lifting surfaces, where the wake detaches and is shed back 
with the reference airspeed, trailing circulation variations without sustaining any pressure load [37]; 
its trajectory then represents a streamline at any time and may generally become part of a nonlinear 
iterative solution where roll-up occurs due to downwash effects, if it is not prescribed a priori (as is 
in fact typical for most practical applications [53,63], especially when characterised by moderately 
unsteady flow). Thus, Morino’s BEM is able to deal with arbitrary 3D geometries in a unified manner, 
reducing the effort of the abstraction process and easing the integration with complex structural models. 
Although the original theory is valid for arbitrary motion in time or frequency domain, only harmonic 
oscillations in the frequency domain are here considered as suitable for flutter analysis. 

The current implementation of Morino’s BEM is described in previous works [64] and is equivalent 
to an appropriate combination of doublets (for lifting bodies and wake) and point sources/sinks 
(for thickness effects and non-lifting bodies). The 3D geometry is approximated with flat quadrilateral 
panels that follow the local wing surface and the aerodynamic potential is assumed constant over 
them, with an analytic expression for their mutual induction. Note that similar flow conditions can be 
modelled by either vortex or doublet distributions and a quadrilateral doublet element is equivalent to 
a vortex ring placed at the panel edges [65]; higher-order methods have also been formulated in order 
to improve accuracy and computational efficiency for complex geometries and configurations, but 
they are more prone to becoming ill-posed and typically not necessary as long as enough quadrilateral 
panels refine the aerodynamic grid for lower-order methods to capture high pressure gradients while 
still staying away from ill-conditioning [66]. Through the wake being assumed as flat, rigid and 
parallel to the free-stream velocity in order to perform rigorous comparisons with Nastran’s DLM and 
Theodorsen’s theory without loss of generality (in the absence of aircraft fuselage and empennages [67]), 
a linear system is obtained: 

T,o* = Ty, (23) 


where frequency-dependent T, and Ty are the matrices of aerodynamic influence coefficients (AIC), 
whereas o* and # are the aerodynamic potential and normal wash associated with the k-th mode 
shape of the structural displacement [68], having mapped the latter onto the aerodynamic grid with 
an effective implementation of the infinite plate spline (IPS) method [69]. The elements Q),; of the 
generalised aerodynamic forces (GAF) matrix Q(s) are then calculated from the work done by the 
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aerodynamic pressure C, due to the k-th mode on the displacement defined by the h-th mode @" for a 
prescribed range of reduced frequencies [40], namely: 


N ‘ 
Qnk = Lop} nj AiC,(o°*), (24) 
i=l 


where n; and A; are the normal vector and area of the i-th aerodynamic panel, respectively. 

The NACA0002 aerofoil [26] is suitably adopted to obtain a baseline 3D representation of Loring’s 
wing surface (see Figure 2); this particular choice of thin symmetric aerofoil is consistent with the 
original work with no loss of generality, as the parametric sensitivity of the aeroelastic stability 
boundary with respect to aerofoil thickness was already found to be poor for small perturbations 
of subsonic potential flow (i.e., in the absence of strong shock waves and/or significant flow 
separation) [70]. 





Figure 2. Aerodynamic panels model of Loring’s wing. 


4.3. Aeroelastic Model 


After obtaining the natural vibration modes, the GAF are computed with Morino’s BEM, and 
the continuation method [50] is used to solve the aeroelastic stability eigenproblem and trace the root 
locus. In order to avoid its costly re-computation during the stability analysis, the GAF matrix is 
approximated by a rational expression in which the non-linear dependency on the complex reduced 
frequency s appears explicitly. Here the matrix fraction approach (MFA) [71] is used and preferred 
to the rational function approximation (RFA) [72], since the former exhibits higher accuracy than the 
latter for the same number of poles [73]; yet, both methods provide the analytical continuation of the 
GAF for complex reduced frequency, increasing the accuracy of damped solutions with respect to the 
p — k method. According to MPA, the GAF matrix is expressed as fraction of matrix polynomials: 


M42 \/M  \7! 
Qr (y Ni (x Ds ; (25) 
i=0 i=0 


and the accuracy of the approximation increases with increasing the number of poles M, which is 
equal to the size of the state-space system [74]; in particular, three poles have been used in this 
proof-of-concept work. The approximation matrices D; and N; are both obtained by solving 
a least-square problem in which the distance from the GAF samples Q(ik) is minimised. 
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The continuation method hereby adopted provides a straightforward and efficient tracking 
technique without using any correlation function, such as the modal assurance criterion (MAC) [75]. 
For this reason, the method is rather insensitive to the number of poles used for the aerodynamic 
finite-state approximation and is able to distinguish between actual and artificial aerodynamic states by 
construction. In the continuation method, the aeroelastic equations are differentiated with respect to the 
free-stream airspeed, thereby resulting in a system of ODEs which are solved with a predictor—corrector 
integration schema starting from an initial airspeed value for which a true solution is known (e.g., 
the natural vibration frequencies for U = 0). 

For validation purposes, flutter analysis is also carried out with Nastran’s DLM directly: a flat 
plate aligned with the free-stream is used as lifting surface and the aerodynamic panels are defined with 
CAERO1 and PAERO!I cards; the distribution and number of such panels follow established guidelines 
and convergence studies (see Appendix B). The mode shapes are mapped onto the aerodynamic grid 
with the IPS method [76] using the SPLINE1 card and the GAF matrix is then generated for sixteen 
reduced frequencies specified in the MKAERO1 card, suitably ranging from k = 0 to k = 1 (with 
logarithmic-like spacing) as per literature studies and common practice [53,63]; a cubic spline is then 
exploited to interpolate the GAF therein. For the dynamic aeroelastic stability analysis, the p — k 
method as available in the FLUTTER card is used: the aerodynamic damping is approximated as 
the imaginary part of the GAF matrix computed for (undamped) harmonic motion, limiting results’ 
accuracy to the case of lightly damped aeroelastic solutions [12]. The p — k approximation of the 
aeroelastic eigenproblem [77] is then solved for all the free-stream speed values defined in the FLFACT 
card, according to the reference length and air density specified in the AERO card. After damping 
ratios are obtained in the given speed range for all modes, the (undamped) flutter point is determined 
by linear interpolation. 


4.4, Aero-Structural Parametric Derivatives 


To calculate the sensitivity of the flutter boundary with respect to any aero-structural parameter, 
the derivatives of the aeroelastic matrices are also necessary in the first place. In particular, 
the derivatives of structural eigenvalues and eigenvectors with respect to structural parameters, such as 
wing material properties (e.g., density or elastic modulus) and element properties (e.g., moments of 
inertia or distance between elastic and inertia axis), are obtained with Nastran SOL200. The design 
variable label and initial value are defined with the DESVAR card, which is connected to the relative 
bulk data entry by DVMREL1 or DVPREL1 cards for material or element properties, respectively; 
derivatives are hence computed for the structural responses defined by DRESP1 cards, which are the 
first n structural eigenvalues and eigenvectors. 

The GAF matrix is differentiated with respect to both reduced frequency and design variables: 
according to MFA, the finite-state approximation allows the analytic differentiation with respect to the 
complex reduced frequency s, whereas the method presented in previous works [73] was used for the 


derivatives with respect to the design variables p. In particular, by defining the functions P (ot, p) 
and R G p) as: 
P= Qn (o%, p), R = Ty ok — Ty, (26) 


a sub-differentiation process is set up, and depending on the number of design variables and mode 
shapes to be considered [11], the derivatives may be obtained by either a direct approach: 


(PP, aPde = ARdo_R - 
dp op oodp’ dodp op’ 
or an adjoint approach: 
dP oP aR dR\* aP\* 
dp ap ap’ (ac) 8=- (55) - " 
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In computing the partial derivatives, it is important to note that the imaginary component 
is introduced only by the normal wash @ and the wake’s influence coefficients inside matrix T;; 
therefore, the partial derivatives of the steady (real) contribution to the GAF are computed via 
complex step [78] and the partial derivatives of the unsteady (complex) contribution are then 
analytically assembled. 


5. Results and Discussion 


Loring’s uniform slender thin wing [25] is hereby considered, as both experimental wind-tunnel 
data and numerical results assuming two-dimensional incompressible potential flow are available 
and all results can hence be explained from both physical and mathematical standpoints; 
moreover, this fundamental benchmark embeds the full complexity of the aeroelastic problem without 
introducing detrimental uncertainties. The wing’s chord is 0.305 m and the semi-span is 2.057 m, 
giving an aspect ratio 13.5. The inertial axis lays at 42.3% of the chord, with mass 8.05 kg/m and mass 
moment of inertia 0.0471 kg-m. The wing root is clamped at the elastic axis, with flexural stiffness 
677.3 N-m? and torsional stiffness 1018.9 N-m? placed at 30% of the chord. 

The coupled natural vibration frequencies of the first bending and torsion modes were observed 
at 1.29 and 18.1 Hz, respectively; that of the second bending mode was detected in between at 7.7 Hz. 
With a fluid density 1.11 kg/m/°, the flutter speed and frequency were experienced at 90.3 m/s and 
10.2 Hz, which give a reduced frequency around 0.11; the subsonic flow may then be considered 
as moderately unsteady. A Euler—Bernoulli beam model calculated the first three coupled vibration 
frequencies as 1.29, 7.65 and 17.98 Hz; once coupled with two-dimensional unsteady aerodynamics 
for a flat plate using standard strip theory, flutter was predicted at 90.7 m/s and 9.2 Hz with good 
accuracy. 

In order to elaborate on the literature results and visualise the flutter mechanism, the aeroelastic 
stability analysis of Loring’s wings is first performed with the same assumptions as in the original 
publication and the approach here described in the Problem Formulation; the p — k method (with 
Theodorsen’s exact function [35]) has consistently been used and results have been cross-verified 
against the state-space representation (with a common two-term RFA of Theodorsen’s function) [21]. 
By employing the first two bending modes and the first torsion mode, the respective coupled vibration 
frequencies are calculated as 1.21, 7.59 and 17.91 Hz; once still coupled with two-dimensional unsteady 
aerodynamics for a flat plate using standard strip theory, flutter is consistently predicted at 91.15 m/s 
and 9.2 Hz (which is indeed an excellent approximation, as confirmed by a modal convergence study 
up to the first three bending and torsion modes). The evolution of the aeroelastic system’s eigenvalues 
is presented in Figure 3 and confirms static divergence to arise well beyond dynamic flutter (i.e., 
Uy >> Uf). 
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Figure 3. Real (left) and imaginary (right) parts of the aeroelastic eigenvalues for Loring’s wing in 
incompressible flow. 
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Due to the remarkable agreement between measurements and simulations, the assumptions of 
slender beam structure and two-dimensional incompressible potential flow in the latter are hence 
deemed justified. Note that the instability mechanism involves second bending and first torsion 
modes directly, but the first bending mode is indirectly essential for their coupling to occur before that 
between first bending and first torsion modes (see Appendix A). 


5.1. Aeroelastic Analyses 


The aeroelastic stability of Loring’s wing is investigated and compared using the proposed 
lower and higher-fidelity models, while still assuming incompressible flow. The analyses encompass 
parametric variations of the aero-structural properties, quantifying their effects on the aeroelastic 
stability boundary and critically assessing the differences in the multifidelity results from both 
theoretical and computational perspectives, for possible practical applications in airplane design 
and optimisation adopting hybrid strategies. 

Figure 4 shows the aeroelastic stability analysis from the higher-fidelity model, focusing on the 
same flutter mechanism as found earlier: the first torsion mode becoming unstable and extracting 
energy from the airflow through the coupled second bending mode, flutter is found at 94.77 m/s and 
10.04 Hz in excellent agreement with the experimental results and exhibits a slightly higher vibration 
frequency than the theoretical one. The aeroelastic stability analysis performed with Nastran and its 
embedded DLM for the flat wing is also presented and confirms the higher-fidelity results based on 
Morino’s BEM with NACA0002 aerofoil, flutter being found at 94.44 m/s and 10.38 Hz. Note that both 
sets of code share the same structural model, and the coupled natural vibration frequencies in the void 
are 1.21, 7.55 and 21.03 Hz for the first three bending modes, whereas it is 17.88 Hz for the first torsion 
mode; these first four modes have been used in all analyses (see Appendix B). As expected from the 
wing aspect ratio being relatively large, three-dimensional downwash effects are actually moderate 
and the (beneficial) steady ones on the (attenuated) airload distribution are roughly compensated by 
the (detrimental) unsteady ones on the (accelerated) airload evolution [21]; thus, the flutter speed is 
just slightly higher than that predicted assuming unsteady two-dimensional flow. 


Re(A) [rad/s] 
Im(A) [rad/s] 
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U [m/s] U [m/s] 





Figure 4. Higher-fidelity real (left) and imaginary (right) parts of the aeroelastic eigenvalues for 
Loring’s wing in incompressible flow. 


Still focusing on the flutter mechanism, Figure 5 then shows the aeroelastic stability analysis from 
the lower-fidelity model: despite its inherent simplicity, the latter predicts flutter quite accurately at 
92.1 m/s and 9.09 Hz, thereby proving the effectiveness of such an idealisation (analogous stability 
diagrams have indeed been obtained by retaining only steady aerodynamic terms in the higher-fidelity 
model, for cross-verification purposes). The coupled natural vibration frequencies in the void are 1.21 
and 7.59 Hz for the two bending modes, and 17.91 Hz for the torsion mode, in excellent agreement 
with the higher-fidelity ones. The aeroelastic response is then metastable until the first instability 
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occurs, as no damping is provided by either the elastic structure or static aerodynamics; to this respect, 
it is worth mentioning that the lift-derivative correction for the steady three-dimensional downwash 
(which tends to increase the flutter speed, by providing lower airload) incidentally compensates the 
lack of aerodynamic lag (which tends to decrease the flutter speed, by neglecting wake inflow). Note 
that, although the flutter reduced-frequency is rather high for steady aerodynamics to be rigorously 
applicable, the latter was still adopted in order to exacerbate the difference between lower and 
higher-fidelity models and stress the multilevel approach for more general conclusions. 


2-nd Bending 
1-st Torsion 


2-nd Bending 
1-st Torsion 


Re(A) [rad/s] 
Im(A) [rad/s] 
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Figure 5. Lower-fidelity real (left) and imaginary (right) parts of the aeroelastic eigenvalues for 
Loring’s wing in incompressible flow. 


5.2. Sensitivity Study 


Following previous studies on the aeroelastic stability boundary of slender cantilever wings 
and its sensitivity [27,28], linear FSI models give nonlinear flutter trends as relevant aero-structural 
parameters are perturbed. The reciprocal positions of aerodynamic, elastic and inertial axes being often 
constrained by the available space inside the wing as well as the chosen structural layout and systems 
arrangement, the flutter point’s sensitivity to the wing’s structural properties is individually explored 
by changing its material density 0m (which scales m and jy) and elastic modulus E (which scales EI 
and GJ), whereas varying the wing’s semi-span | (which modifies the lift coefficient derivative Cj /,) 
alters its geometry and affects both structural and aerodynamic properties at the same time. 

From Figure 6, it can be seen that changing the wing’s material density alters all natural frequencies 
through the inertial properties and hence the flutter frequency, but has marginal/no influence on the 
flutter speed; all symbols give the actual nonlinear percentage variation of the flutter point, whereas all 
lines with corresponding colour draw its normalised linear prediction. In particular, a 1% increment in 
the material density causes about a 0.5% decrement in the flutter frequency; the negligible variation of 
the flutter speed is confirmed by previous works [27] as being due to the (small) beneficial effect of an 
increase in sectional mass being compensated by an (almost) equal and opposite detrimental effect of 
an increase in mass moment of inertia. This outcome may indeed help while minimising the aeroplane 
weight, as reducing the material density does not significantly affect the present flutter boundary due 
to wing bending-torsion instability (especially in the absence of lumped masses). However, it shall be 
recalled that varying the wing inertia might cause other types of aero-structural instabilities to arise at 
the aircraft level due to modal coalescence, and monitoring the variation of the flutter frequency is then 
important to preventing potential resonance. It is also worth mentioning that the small higher-fidelity 
effect on the flutter speed is due to the unsteady airload being frequency dependent through the 
lift-deficiency function and the related aerodynamic lag. 
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Figure 6. Flutter speed (left) and frequency (right) parametric variation and sensitivity to changes in 

the wing’s material density. 


As per Figure 7, changing the wing’s material elastic modulus alters all natural frequencies 
through the stiffness properties and hence the flutter speed and frequency; in particular, a 1% increment 
in the elastic modulus induces about a 0.5% increment in all the latter, leaving the flutter reduced 
frequency practically unchanged. This outcome is quantitatively confirmed by previous works [27] 
as being due to the detrimental effect of an increase in flexural stiffness being much lower than 
the beneficial effect of an increase in torsional stiffness. The striking agreement between lower and 
higher-fidelity results is mostly driven by the respective structural models being equivalent (as seen 
from the close agreement between the natural vibration frequencies), while the different aerodynamics 
play a minor role (as the airflow is moderately unsteady). The quasi-linear trend of the percentage 
variations reveals a rather large confidence interval for the normalised sensitivity; however, the 
dimensional counterpart of the latter to be used by optimisation routines follows a highly nonlinear 
pattern (note that the explicit lower-fidelity expression given in Appendix A for the torsional static 
divergence speed provides a straightforward theoretical check). 
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Figure 7. Flutter speed (left) and frequency (right) parametric variation and sensitivity to changes in 
the wing’s material elastic modulus. 


Finally, Figure 8 shows that changing the wing semispan has a large effect on the stability 
boundary through the stiffness properties as well as the aerodynamic loads and hence on the flutter 
speed and frequency. In particular, a 1% increment in the wing semispan induces about a 0.8% 
decrement in the critical speeds and about a 1.6% decrement in the flutter frequency; note that 
analogous trends for the variations and orders of magnitude for the sensitivities were obtained in 
previous studies on similar slender wings [28], as a qualitative means of validation. Lower and higher- 
fidelity results are still in very good agreement and differences are mainly due to unsteady downwash 
effect, as the aspect ratio of the wing changes considerably with the span and so does the sectional airload. 
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Figure 8. Flutter speed (left) and frequency (right) parametric variation and sensitivity to changes in 
the wing semispan. 


For the sake of completeness, it is worth mentioning that variations and sensitivities of the 
divergence speed (see Appendix A) with respect to the same aero-structural parameters considered 
above were consistently found to have trends and orders of magnitude very close to those pertaining 
the flutter speed. In particular, it is observed that changing the wing’s geometry has a much larger 
impact on the flutter boundary than changing the material properties (especially density); this is 
particularly true for the flutter frequency, the variations of which exhibit a more significant nonlinear 
behaviour and local curvature. 

Due to the sufficient mutual separation of the natural vibration modes, it is also worth stressing 
that the instability mechanism did not change across the parametric variations, and the derivatives of 
the flutter stability boundary have accurately been calculated by the lower-fidelity model at a (marginal) 
fraction of the computational complexity and costs. In particular, semi-analytical solutions drastically 
reduced the latter, being almost instantaneous and demanding minimal pre- and post-processing (if any 
at all) while granting an enhanced theoretical understanding. Thus, the healthy combination of lower 
and higher-fidelity models enables efficient multidisciplinary exploration of a large design variable 
space for innovative aeroplane concepts and configurations; further numerical savings may still be 
obtained by exploiting reliable surrogate models for the higher-fidelity solutions, with an effective 
synthesis of the underlying complexity [20]. 


6. Conclusions 


Within the context of aircraft multidisciplinary design and optimisation, a comparative sensitivity 
study for the bending-torsion flutter instability of flexible aircraft wings in subsonic flow has been 
presented. Analytical models and numerical tools with different complexities and fidelities have 
been used, in view of possible practical applications exploiting multilevel approaches within the 
conceptual and preliminary MDO phases. Parametric studies have been performed where the effects 
of varying the wing’s aero-structural properties on the aeroelastic stability boundary have been 
quantified and critically assessed from both theoretical and computational perspectives. When the 
natural vibration modes of the wing are well separated from all other natural vibration modes of 
the aircraft and the aeroelastic instability mechanism does not change in nature, an efficient hybrid 
strategy is then recommended where the flutter analysis is performed using higher-fidelity approaches, 
whereas the sensitivity analysis of the flutter boundary is performed using lower-fidelity approaches, 
thereby improving theoretical understanding and reducing computational costs while retaining 
accuracy. Future works are encouraged to investigate additional effects (e.g., wing sweep and flow 
compressibility) and increased complexity (e.g., presence of a control surface) or perform the full MDO 
of a flexible aircraft wing, effectively exploiting the proposed multifidelity strategy. 


Aerospace 2020, 7, 161 15 of 22 


Author Contributions: M.B. derived the analytical model and results; F.T. performed the numerical simulations; 
the authors then wrote the respective parts of the manuscript. All authors have read and agreed to the published 
version of the manuscript. 


Funding: This research received no external funding. 


Acknowledgments: The authors would like to thank Ranjan Banerjee at City University of London and Rakesh 
Kapania at Virginia Polytechnic Institute and State University for the precious insights into their previous works. 


Conflicts of Interest: The authors declare no conflict of interest. 


Abbreviations 

Symbols 

A aerodynamic panel area 

C section chord 

C section lift 

Cia section lift derivative 

Cry wing lift derivative 

Cp pressure coefficient 

Cr Theodorsen’s function 

[C] generalised damping matrix 

D aerodynamic approximation matrix (fraction denominator) 
e semiperimeter-to-span ratio 

E section Young’s elastic modulus 

f cross-projection of first bending and first torsion modes 
{F} generalised load vector 

g cross-projection of second bending and first torsion modes 
G section shear elastic modulus 

h section flexural (plunge) displacement 

jake Hankel’s functions of the second type and n-th order 

I section flexural area moments of inertia 

J section torsional area moments of inertia 

k reduced frequency 

k equivalent spring stiffness 

[K] generalised stiffness matrix 

l wing semi-span 

AL section aerodynamic force 

m section mass 

AM section aerodynamic moment 

|M] generalised mass matrix 

n aerodynamic panel normal vector 

N aerodynamic approximation matrix (fraction numerator) 
p design parameter 

Q generalised aerodynamic forces matrix 

r squared ratio of second and first flexural vibration frequencies 
s complex reduced frequency 

[S] system matrix 

t time 

T aerodynamic influence coefficients matrix 

{u} eigenvector 

U horizontal airspeed 

V vertical airspeed 

x chordwise coordinate 

y spanwise coordinate 

Ww section vertical displacement 
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Greek 
ni 


{x} 


es ae SS oS A es He 


Subscripts 


SDN FAW LY 


Acronyms 


AC 
AIC 
BEM 
CFD 
CG 
CP 
CSRD 
DLM 
EA 
FEM 
FSI 
GAF 
IPS 
MAC 
MC 
MDO 
MFA 
MST 
ODE 
PDE 
OST 
RFA 
ROM 
SST 
TST 


angle of attack 

generalised coordinates 

aerofoil thickness ratio 

natural vibration mode shape 
flexural natural vibration constant 
wing aspect ratio 

aerodynamic load scaling function 
eigenvalue 

section mass moment of inertia 
torsional natural vibration constant 
section torsional (pitch) displacement 
normal wash matrix 

air density 

material density 

aerodynamic potential matrix 
reduced time 

natural vibration frequency 


aerodynamic 
critical 
flutter 
divergence 
flexural 
structural 
torsional 


aerodynamic centre 
aerodynamic influence coefficient 
boundary element method 
computational fluid dynamics 
centre of gravity 

control point 

closely-spaced rigid diaphragm 
doublet lattice method 

elastic axis 

finite element method 
fluid-structure interaction 
generalised aerodynamic forces 
infinite plate spline 

modal assurance criterion 
mid-chord 

multidisciplinary design and optimisation 
matrix fraction approach 
modified strip theory 

ordinary differential equation 
partial differential equation 
quasi-steady theory 

rational function approximation 
reduced order model 

standard strip theory 

tuned strip theory 
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Appendix A. Aeroelastic Stability of the Typical Section with Steady Aerodynamics 


Providing full control on the results with respect to the specific assumptions of the methods and 
tools employed for the analysis, theoretical formulations grant a clear and complete overview of the 
problem which is essential for any engineering application. Although inherently limited in their general 
capabilities, analytical solutions then offer a wealth of qualitative information and quantitative details 
as well as fundamental insights and rigorous validation for both educational and practical purposes. 

Considering a single mode for the wing bending and torsion, respectively, the aeroelastic equations 
for the Typical Section with steady aerodynamics read [79]: 


m (ht —xcc fb) +k,h = EU*cCy inf, (Al) 


ud —MXxcG (h — xcc fé) +kg@ = EU xaccCr jab, (A2) 


where f = 1 for two-dimensional aerofoils [12], with ¢, = 1 and ¢g = 1. By neglecting the aerodynamic 
load, the coupled bending and torsion natural vibration frequencies are explicitly obtained as: 


1 
Fyws — Fow* + Fy = 0, Wy = oe (F + 4/ FR — Fifa), (A3) 
4 


_ Hw; Ws 
wt mx2.(1— fe) 


where the uncoupled counterparts are resumed whenever inertial and elastic axes coincide (i.e., 


Ht 2 mM 9 2 
Fy = 1, R= om é  — joe+ii+-—x wrt, E A4 
= mee (et (1+ Fee) a] = 


with xccg = 0). Otherwise, the characteristic equation for the flutter determinant provides with the 
metastable boundary and gives the flutter frequency as: 


1 
Ae. He _ [p2 _ 
Paco Pow + Pp =0, wr= DP, (?, + 4/ Ps 4P,P) , (A5) 


and setting the inner radical discriminant equal to zero then gives the flutter speed as: 


1 
D4U; = D2UF + Do = 0, Ur = Dh (-p;  4/ Ds = 4D,Dp)), (A6) 


where the sign before the inner radical shall provide with the smallest positive real value; 
finally, note that setting both w+ and Po equal to zero provides with the torsional static divergence 
speed. In particular, when the aerodynamic centre is ahead of the elastic axis (i.e., with x4c < 0), 


U, = ne. (A7) 
OX ACCC Ly 


regardless the bending stiffness; when the elastic axis is ahead of the inertial axis (i.e., with xcc¢ > 0), 


the latter explicitly reads: 


the flutter condition is given as: 


~ J+ (/pe : — {22% 
Uf = or ( Ds 4D4Dpo D2), Wf = DP,’ (A8) 


where all coefficients are analytical functions of the aero-structural parameters, namely: 


i= kp (ke + FUP?xaccCr/x) ; (AY) 
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Pp=m ke + eu (xc — xcof?) Cra + ky, (1 + mx) ; (A10) 
Py =m |p t+mxec (1-f?) |, (A11) 
Do = | kg + ky (1+ mx2c)| — 4m |p + mxeg (1— f?)| knke, (A12) 


Dz = 2m { (xc — xcof?) kam + ky (1 at mxeg)) | ae, f + Mxec (1 — f°) xacks } CCr/y, (A13) 


y= li (xc - xcof?) Ciya] (A14) 


The formulas above were implemented and compared with literature results for the stability 
boundary of a Typical Section assuming standard strip theory (SST) with steady aerodynamics [12,79]: 
exact agreement was always found and provided rigorous validation, also with respect to the relative 
root locus. The same formulas have then been used to investigate the stability boundary of Loring’s 
wing Typical Section when either the first or the second bending mode is coupled with the first torsion 
mode, employing tuned strip theory (TST) with steady aerodynamics [21]. Without accounting for the 
modal cross-projection (i.e., f = 1 fora “pitch & plunge” apparatus), flutter is calculated at 106.5 m/s 
and 4.32 Hz in the former case whereas at 73.9 m/s and 11.28 Hz in the latter case; when the modal 
cross-projection is considered (i.e., f < 1 for a slender beam), flutter is calculated at 109.7 m/s and 
4.28 Hz in the former case whereas at 139.2 m/s and 9.60 Hz in the latter case. The static divergence 
speed pertains the torsion mode only and is found at 210.2 m/s, regardless the bending modes and 
their cross-projections. 

A direct comparison with the low-fidelity results hence reveals that the interaction between 
bending and torsion modes as well as their cross-projections are essential to reproduce the correct 
behaviour of the flutter mechanism, where the coalescence between first bending mode and torsion 
mode drives that between second bending mode and the latter (which have closer natural vibration 
frequencies) to occur earlier. 


Appendix B. Higher-Fidelity Model Results Convergence Study 


The number of elements for the structural and aerodynamic models was defined according to 
rigorous convergence studies. Following Nastran’s best practice [80] with kingy = 1 the maximum 
reduced frequency employed in the aeroelastic analysis, fifteen DLM panels have uniformly been 
distributed along the wing chord. When adopting Morino’s BEM, thirty aerodynamic panels were 
symmetrically placed along both upper and lower aerofoil surfaces according to the convergence 
study in Figure A1 (left), with a suitable refinement around the leading edge in order to capture high 
pressure gradients. The number of aerodynamic panels in the span-wise direction was determined 
according to another convergence study shown in Figure Al (right) and all higher-fidelity results 
presented in this work were obtained with fifteen panels strips uniformly distributed along the wing 
span, since sufficient to grant a relative error below 1% for both flutter speed and flutter frequency. 
Following convergence studies in previous works [37,64], the wake extends for a hundred chords 
behind the wing and is modelled with a hundred rows of shed panels, the length of which is cubically 
increased with increasing their distance from the trailing edge. Figure A2 (left) shows the indicial 
lift-deficiency function ®(T) (which represents the equivalent of the complex lift-deficiency function 
Cr in the reduced-time domain T) from a unit step in angle of attack, whereas Figure A2 (right) 
depicts the normalised distribution of the steady circulation «(y) (which represents the ratio between 
the sectional circulation in three- and two-dimensional flow within the framework of modified strip 
theory, MST [13]); note that the unsteady lift development approaches quasi-steady theory (QST) 
asymptotically. These results also justify the assumption of two-dimensional flow for the lower-fidelity 
model, within the framework of TST [21] (to this respect, it is worth mentioning the higher-fidelity 
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lift coefficient derivative of the wing is obtained as Cj, = 5.25, in excellent agreement with the 
lower-fidelity analytical estimation Cy /, = 5.21). 


% Difference 








Nchord 





Figure A1. Flutter speed and frequency varying the number of aerodynamic panels along the wing 
chord (left) and span (right). 











T 


Figure A2. Indicial lift-evolution function from a step-change in the angle of attack (left), spanwise 
lift-decay function (right) for Loring’s wing in incompressible flow. 


Figure A3 (left) then shows the evolution of the natural vibration frequencies percentage error with 
varying the number of FEM nodes of the structural beam; the value obtained with the most dense grid 
is used as reference to compute the plotted error. The first four natural vibration frequencies exhibit 
a good convergence behaviour, with the third one (i.e., the first torsional mode) almost insensitive 
to the number of nodes. All higher-fidelity results presented in this work were hence obtained with 
20 nodes uniformly distributed along the wing span, since sufficient to grant a relative error below 
2%. Finally, Figure A3 (right) shows that employing the first four natural vibration modes granted 
convergence of both flutter speed and flutter frequency. 
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Figure A3. Natural vibration frequencies varying the number of FEM nodes along the wing span (left), 
flutter speed and frequency varying the number of natural vibration modes employed in the aeroelastic 
stability analysis (right). 
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Abstract: The German Cluster of Excellence SE*A (Sustainable and Energy Efficient Aviation) is 
established in order to investigate the influence of game-changing technologies on the energy efficiency 
of future transport aircraft. In this paper, the preliminary investigation of the four game-changing 
technologies active flow control, active load alleviation, boundary layer ingestion, and novel materials 
and structure concepts on the performance of a long-range Blended Wing Body (BWB) aircraft is 
presented. The BWB that was equipped with the mentioned technologies was designed and optimized 
using the multi-fidelity aircraft design code SUAVE with a connection to the Computational Fluid 
Dynamics (CFD) code SU2. The conceptual design of the BWB aircraft is performed within the 
SUAVE framework, where the influence of the new technologies is investigated. In the second step, 
the initially designed BWB aircraft is improved by an aerodynamic shape optimization while using 
the SU2 CFD code. In the third step, the performance of the optimized aircraft is evaluated again 
using the SUAVE code. The results showed more than 60% reduction in the aircraft fuel burn when 
compared to the Boeing 777. 


Keywords: aircraft design; aerodynamics; aerodynamic design optimization; blended wing body; 
new airframe technologies 


1. Introduction 


A series of research activities have been carried out all over the world in the recent years to meet 
the challenges of ambitious reduction in CO2, NOx, and noise emission set by aviation authorities, 
such as in Flightpath 2050 [1]. Among which one notable effort is the German Cluster of Excellence 
SE*A (Sustainable and Energy Efficient Aviation). The SE*A program is based on the previous joint 
research project “Energy System Transformation in Aviation, EWL [2]” that has been initiated between 
2016 and 2018 in Germany to identify and investigate possible unconventional energy systems that 
can be used for civil transport aircraft in combination with game-changing aircraft configurations and 
airframe technologies. The SE*A project aims at investigating the influence of new technologies as 
well as new operational scenarios on the sustainability of future transport aircraft. It has been inspired 
by a lot of relevant work that has been carried in recent decades; some are mentioned in the following. 

A blended wing body (BWB) concept for long-range commercial aircraft could lead to a fuel saving 
of 27% as compared to a conventional A380-like tube-and-wing (TAW) configuration, according to the 
research of Boeing [3]. Since the first appearance of the BWB as a potential future commercial aircraft 
configuration, a lot of research related to the BWB has been done. Okonkwo and Smith [4] provided 
a summary of achievements in the design of the BWB and its benefits as compared to conventional 


Aerospace 2020, 7, 87; doi:10.3390/aerospace7070087 www.mdpi.com/journal/aerospace 


Aerospace 2020, 7, 87 2 of 27 


aircraft. Xu and Kroo [5] investigated the benefits of load alleviation and natural laminar flow and 
they concluded that the combination of these two technologies to a Boeing 737-800 aircraft could bring 
a fuel saving of 18%. NASA-MIT D8 “DoubleBubble” concept with boundary layer ingestion (BLI) 
and active load alleviation has conducted a fuel burn reduction of 70.87% as compared to a B737-800 
baseline concept [6,7]. NASA Hybrid Wing Body (HWB) concept (with N + 3 airframe technology 
packages, such as BLI and distributed propulsion system) has conducted a fuel burn reduction of 54% 
as compared to a 777-200LR baseline [6,7]. The “Advanced Truss-Braced Wing” concept proposed by 
NASA and MIT with hybrid-electric propulsion has a 70% fuel burn reduction at a very low range 
condition [8]. Saeed et al. from Cambridge University have designed a flying wing concept with 
laminar flow control and concluded that with 84% of the total wetted area being laminarised, they 
have achieved a 70% fuel savings when neglecting the system penalties [9]. 

Within the EWL project, several aircraft have been initially designed to represent technology 
integrator for these new technologies. In particular, the following technologies have been selected via 
system-level studies [2]. 


1. Active Flow Control: the function of this technology is to actively suck the air from the aircraft 
outer surface to delay the transition of the boundary layer. The skin of each aircraft component is 
split into two segments: a porous sheet and an inner sheet that supports the outer sheet. The inner 
sheet has orifices that suck the air from the boundary layer and delay transition. In each chamber, 
an individual pressure is adjusted by the throttle orifices, so that the pressure difference between 
the outside and the chamber delivers the locally desired amount of mass flow through the surface. 
Figure 1 shows a schematic image of the skin layout and the AFC system for a wing section. 
The applied technology in this project is based on [10,11], which describes numerical approaches 
with active laminar flow control and also describes current progress in this technology. 
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Figure 1. Schematic views of the active suction system [10]. 


2. Active Load Alleviation: the wing active load alleviation consists of several technologies that 
actively reduce bending moment experienced by the wing, so the limit load factor for the aircraft 
can be reduced and, consequently, the wing weight is reduced. Previously, researchers have 
approached the design of active load alleviation systems in different ways. For example, [12-14] 
looked at the alleviation technology from the aircraft design and MDO perspective; [15] looked 
at active load alleviation from the control system perspective; and, [16] demonstrated the 
experimental results of high aspect-ratio wing wind-tunnel test with active load alleviation while 
using piezoelectric control. 

3. Boundary Layer Ingestion (BLD: in this technology, the jet engines are integrated with the aircraft 
body in such a way that the boundary layer of the aircraft body is ingested into the engines, 
which improves Specific Fuel Consumption (SFC). Extensive information about the current state 
of the BLI modeling, analysis, and its benefits from the aircraft design standpoint are provided 
in [17-21]. In addition to BLI, ultra-high bypass ratio turbofan engines that further improve fuel 
consumption [22] are assumed in this project. 

4. |New Materials and Structure Concepts: novel structural concepts and materials are developed 
to improve the aircraft structure in terms of stiffness and weight. Bishara et.al. [23] describes 
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advanced structural design with the integration of active flow control, which is particularly 
important for the presented research. Besides, aeroelastic tailoring is also considered in this project. 


The focus of the current manuscript is to present the multi-fidelity design and analysis of one 
of the reference aircraft, i.e., a long-range BWB aircraft, which is used to demonstrate the benefits 
of combining the above-mentioned technologies for improving the fuel consumption of transport 
aircraft. The studies performed in [2] determined the initial configuration of the aircraft and estimated 
geometric properties that are based on mission requirements and several trade studies, such as the 
wing aspect ratio, taper, sweep, and thickness. This study focuses on aircraft mission analysis and 
optimization while using a multi-fidelity approach. In this manuscript, first, the methods and tools that 
are used to carry out the overall aircraft design and optimization are briefly described. Subsequently, a 
more complete description of analysis methods and the outcome of the initial (conceptual) design of 
the BWB aircraft are presented. Finally, an aerodynamic shape optimization method and its results 
are presented. 


2. Methods and Tools 


As mentioned before, the goal of this research is to investigate the influence of the new technologies 
mentioned above on the fuel consumption of a long-range BWB. To achieve this goal, a multi-fidelity 
design optimization is performed in order to design a long-range BWB aircraft equipped with the 
mentioned technologies. Three open-source tools were used in this research. The overall assessment 
of the aircraft was executed while using the Stanford University Aerospace Vehicle Environment 
(SUAVE) [24]. Aerodynamic analysis and optimization of the aircraft were executed while using the 
SU2 Computational Fluid Dynamics (CFD) code [25]. The SUAVE code was connected to the OpenVSP 
code [26] for automatic CAD modeling and CFD-mesh generation for SU2 from the aircraft geometry 
defined in SUAVE. 


3. Conceptual Design and Assessment of the BWB Aircraft 


3.1. Top Level Aircraft Requirements and Initial Design 
Table 1 lists the top-level aircraft requirements of the reference long-range aircraft SE* A-LR, which 


were derived from the same category transport aircraft such as Boeing 777/787 or Airbus A330/350. 


Table 1. Top-level aircraft requirements of reference long-range aircraft Sustainable and Energy Efficient 
Aviation-LR (SE*A-LR). 


Parameter Unit Value 
Design range NM 8099 
Design passenger number - 300 
Cruise Mach number - 0.85 
Maximum Mach number - 0.92 


Cruise altitude m 10,600 
Take-off field length (TOFL) m <2200 
Landing distance m <1966 


The initial design of SE*A-LR is based on the EWL project [2], where the thrust to weight ratio 
and wing loading values were determined by top-level aircraft requirements and design specifications 
derived from existing long-range transport aircraft, such as Boeing 777. The outer wing of the BWB 
SE*A-LR was designed via optimizing planform parameters, such as wing aspect ratio, leading-edge 
sweep, thickness to chord ratio, and taper ratio for the objective function of MTOW for a flight range 
of 15,000 km while using the NOMAD optimizer [27]. The center body of SE*A-LR was designed to 
carry 300 passengers using a multi-bubble concept. The supercritical DLR F15 airfoils were used for 
the initial design with different thickness to chord ratios at the center body and outer wing sections. 
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The three-dimensional (3D) view of the initial design and comparison of geometric planforms between 
the SE7A-LR and B777 and their wing properties are shown in Figures 2 and 3, respectively. 





Figure 2. Three-dimensional (3D) view of the initial reference long-range blended wing body 
(BWB) aircraft. 





EWL-LR reference concept 
Wing area: 363 m? 

Wing AR: 9.2 

Wing span: 57.8m 

Wing taper: 0.25 

Wing leading edge sweep: 29 





Figure 3. Geometric comparison between long range aircraft SE*A-LR and B777 outer wings 
planforms [2]. 


It has to be noted that, for the initial design of SE7A-LR, hybrid laminar flow control (80% of the 
blended wing area is laminarised), boundary layer ingestion (5% improvement in propulsive efficiency 
due to boundary layer ingestion), advanced structure design with composite materials (20% structure 
weight reduction as compared to baseline), as well as active load alleviation technologies (load factor is 
reduced to 1.5 g) were included [2]. Table 2 gives a summary of the new airframe technologies applied 
to the reference long-range aircraft SE*A-LR [10,12,22,23]. 


Table 2. A summary of new airframe technologies applied to reference long-range aircraft SE? A-LR. 


SE*A Technology Assumption 
Active Laminar Flow Control 80% of the aircraft area is laminarized 
Advanced structures 20% structure weight reduction 
Active Load Alleviation Limit load factor of 1.5 
Boundary layer ingestion 5% improvement in propulsive efficiency 


Lateral—Directional stability and control are carried by two winglets that feature rudders. The tails 
were sized based on static stability requirements presented in [28]. Based on [28], and are required for 
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sufficient stability of the BWB. The aerodynamic derivatives required for static stability analysis and 
tail sizing were computed using the Vortex Lattice code AVL [29]. For the SE7A-LR, and are computed 
equal to and respectively at the cruise conditions. This satisfies sufficient stability requirements. Table 3 
provides a summary of the geometric parameters of the SE*A-LR wing. 


Table 3. Wing and vertical tail geometric properties. 





Geometric Parameter Wing Winglet 
Span (m) 57.80 3.00 
AR (m) 7.40 2.00 
Taper ratio (m) 0.065 1.00 
Root chord (m) 30.60 1.5 
Quarter chord outer wing sweep 20.00 30.00 
(deg) 


SE*A-LR features four engines located on top of the top of the inboard wing portion close to the 
wing trailing edge, as shown in Figure 2. Safran CFM 56-7B22 engines were selected for the aircraft 
based on the required thrust. Table 4 shows the specifications of CFM 56-7B22. 


Table 4. Safran CFM 56-7B22 engine specifications [30]. 


Engine Specifications CFM 56-7B22 
Length (cm) 250.8 
Width (cm) 211.8 
Height (cm) 182.9 
Dry weight (kg) 2386.0 
Take-off thrust (kN) 91600 
Take-off Thrust-specific fuel consumption (g/kN/s) 10.1 


3.2. Assessment of the Benchmark Aircraft and the Initial Design Using SUAVE 


Because the initial design of the BWB aircraft in the EWL project has been done with an in-house 
tool (different from SUAVE), in the first step, we assessed the initial configuration of the SE2A-LR, 
including the game-changing technologies in SUAVE, which is an open-source, object-oriented aircraft 
design environment programmed in Python language with good flexibility, composability, and 
extensibility [24,31]. It enables multi-fidelity analyses of arbitrary aircraft and propulsion systems (both 
conventional and unconventional aircraft concepts as well as propulsion systems). The performance of 
desired components in SUAVE is calculated using individual design or analysis modules with multiple 
fidelities for different cases. For aerodynamics, both build-up methods (including AVL Vortex Lattice 
code for induced drag calculation), and higher-fidelity CFD approach (SU2 as solver and OpenVSP for 
CFD meshes) are used. The interfaces with Gmsh for generating a mesh for SU2 are also available 
in SUAVE [32]. Currently, empirical and statistical methods or surrogates are used for structure and 
weight estimation. A modular “energy network” has been implemented in SUAVE based on analytical 
methods that are used for both gas-turbine and electric energy systems (electric motor, fuel cells, 
batteries, etc.). The aircraft mission in SUAVE is analyzed by iteratively solving the equations of motion 
with a segment-based architecture [24]. By comparing SUAVE analysis results for Boeing 737-800, 
Embraer E-190, Concorde, and Boeing SUGAR Ray BWB with literature, the SUAVE tool has shown 
good accuracy for a wide range of transport aircraft [24]. The aircraft geometry in SUAVE is described 
using representative parameters that can be used for simple aerodynamic/structural analyses, such 
as VLM and Beam Theory. By using additional geometry converter, such as OpenVSP, the aircraft 
seometry parameterization in SUAVE can be further used to generate CFD meshes for high fidelity 
aerodynamic studies [33]. 
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Initial assessment using SUAVE was performed using a low/medium fidelity model, where 
AVL was used for induced drag, and semi-empirical formulations were applied to estimate parasitic, 
compressibility, and miscellaneous drag. The detailed Operating Empty Weight (QWE) mass breakdown 
was also estimated while using semi-empirical formulations within SUAVE. The results were comparable 
with those of the EWL project [2]. Appendix A summarizes SUAVE analysis for B777 and SUAVE 
simulation results for the SE7A-LR concept. Table 5 shows the comparison between SE7A-LR and 
baseline B777, a similar mission for which was also simulated using SUAVE. From this table, one 
can observe that the new BWB aircraft MTOW was reduced by 58.5%, and block fuel was reduced 
by 73.6% as compared to the Boeing 777. Aerodynamics performance, especially the L/D, has been 
significantly improved due to the BWB configuration and laminar flow control technology. During 
the mission, L/D reaches the value of 48.0, while the B777 model reaches L/D of 21. If the SE7A-LR 
concept is compared against existing BWB concepts summarized in [34], significant improvements in 
aerodynamic efficiency are also observed. As predicted by [34,35], the range of maximum L/D does 
not exceed 28 for existing BWB concepts, see Figure 4. Therefore, the improvement of aerodynamic 
efficiency of the SE*A-LR when compared to the most efficient BWB is equal to 71.4%. Such a difference 
between the aerodynamic efficiency of two BWB aircraft is the result of an extended laminar flow 
up to 80% of the chord. In addition to that, the fuel efficiency of the SE*A-LR was compared to the 
BWB aircraft that was presented in [36], which has been designed for similar top-level requirements. 
However, the BWB of [36] was designed without considering the novel control technologies of SE*A. 
Due to this fact, the difference in fuel efficiency between the BWB in [36] and the SE*A-LR was equal 
to 73.6%. 


Table 5. Comparison of key aircraft parameters of reference long-range aircraft SE*A-LR and B777. 


2 ° 
Parameter as ae SE2A-LR (SUAVE) _B777 greeks 
MTOW (kg) 132,268.0 144,308.0 347,452.0 —58.5 
OWE (kg) 78,472.0 82,484.0 145,150.0 —44.2 
L/Dmax 45.5 48.0 21.0 127.0 
Cruise average L/D 38.33 34.0 18.5 83.8 
Block fuel (kg) 24,623.0 28,824.0 109,290.0 —73.6 
Sea level static thrust (KN) 400.0 442.3 1026.0 —56.9 
Fuel Efficiency (kg/seat/100 km) 0.55 0.64 Pol —76.5 
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Figure 4. Summary of maximum L/D for various types of aircraft [32]. 


Appendix A provides details of the mission analysis of SE*A-LR, including the information of 
altitude and weight changes through flight mission as well as the velocity data along the flight mission. 
It has been noticed that there are some deviations between the results of the initial analysis of SE*A-LR 
and the analysis using SUAVE, as shown in Table 5. Such a difference appeared due to a more detailed 
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mission analysis approach in SUAVE when compared to the in-house tool, which presents a more 
accurate estimate of potential fuel savings of the BWB. 


4. High Fidelity Aerodynamic Analysis 


A higher fidelity aerodynamic analysis and optimization was performed after the initial analysis 
of SUAVE using low to medium fidelity methods. For high fidelity aerodynamic analysis, a Stanford 
University Unstructured (SU2) CFD tool [25] was used. In SU2, the finite volume method is 
employed to discretize the Euler and RANS equations, with both explicit and implicit methods 
available for time integration. Using the Free-Form Deformation (FFD) technique, SU2 computes the 
deformation of two-dimensional (2D) and 3D geometries within the computational mesh. Besides, the 
adjoint implementation of the Euler and the Navier-Stokes equations enables efficient gradient-based 
aerodynamic shape optimization using SU2. More details on SU2 can be found in the literature [25,37]. 

The NASA OpenVSP tool [26] is used to link SUAVE and SU2. In OpenVSP, the aircraft geometry 
is described in XML format, which can be easily connected to high fidelity analysis tools, such as SU2. 
For example, the surface triangulations in OpenVSP can be read by the open-source Gmsh tool with 
MSH output format. More recently, OpenVSP has also enabled the capability of creating CFD and 
FEM mesh directly from geometry data (.VSP3 file). Additional information on OpenVSP is available 
in references [12]. Via an automatic link between SUAVE and SU2, as illustrated in Figure 5, the 
aerodynamic performance of the aircraft was evaluated using Euler analysis in SU2 (induced and wave 
drag) in order to improve the accuracy of the analysis. The low-fidelity flat plate analogy [38] method 
was used for the viscous drag estimation. 


Generation of geometry Visualization of geometry Mesh generation 


OpenVSP 





Analysis 


Figure 5. Sequential build-up from Stanford University Aerospace Vehicle Environment (SUAVE) 
geometry and to SU2 CFD analysis. 


To extract the compressibility drag from Euler analysis, the following formulation was used 


CDeomp = CD —- Colm =03 (1) 


where Cp,,,,,, is the drag component due to compressibility and Cpljy = 03 1s the drag at M = 0.3, the 
Mach number where compressibility effects are negligible. An additional SUAVE script was written to 
extract compressibility drag from SU2 and include it in the analysis as a surrogate. For the surrogate 
model, instead of the default Gaussian process surrogate model in SUAVE, the Surrogate Modeling 
Toolbox (SMT) developed by the MDO lab of the University of Michigan [39] was integrated in SUAVE. 
The RMTB B-spline method was used to fit the data. The SMT-toolbox was chosen because of its 
superior fitting capabilities when compared to the default Gaussian process method with SUAVE. 
A comparison of the lift and compressibility drag coefficients between the two methods are shown 
in Figures 6 and 7. From these figures, one can observe that the RMTB method presents much better 
flexibility to accurately fit both lift and drag compared to the Gaussian process. 
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Figure 6. Comparison between the RMTB and Gaussian process surrogates for the compressibility 
drag coefficient. 
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Figure 7. Comparison between the RMTB and Gaussian process surrogates for the compressibility 
drag coefficient. 


The aerodynamic results that were calculated with SU2 were used to build up surrogates for 
SUAVE overall aircraft aerodynamic performance estimation. Figure 8 shows the airframe geometry 
and an example CFD mesh of half of the SE7A-LR aircraft geometry. It is important to note that the 
effect of winglets is excluded in the present CFD analysis. Such a decision was made due to limitations 
in the definition of the FFD-box (see Section 5) and to adequately compare planforms before and after 
the aerodynamic optimization. It is assumed that the winglet aerodynamic problem will be directly 
solved during the next design iteration and the wing optimal planform will be updated again. 





Figure 8. Airframe geometry of SE*A-LR aircraft visualized in OpenVSP. 
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A mesh convergence analysis has been performed to accurately capture the compressibility effect. 
The results of the drag coefficient sensitivity to mesh resolution at zero angle-of-attack and M = 0.85 
are shown in Table 6. Here, an extrapolated value was estimated using the following procedure. 


Table 6. Mesh convergence study for the cruise configuration in SU2 (« = 0 deg, M = 0.85). 


Mesh Number Cell Cp Error (%) 
1 343,612 0.01733 9.02 
2 660,528 0.01673 5.25 
3 835,000 0.01670 5.03 
4 1,012,683 0.01665 4.72 
Extrapolation CO 0.01590 0.00 


The computational domain volume is assumed as 
=3 
V = Nh (2) 


where N is the number of cells and h is the average cell side length. Subsequently, the surface cell area 


is described as 
) V 1 


=e ee (3) 
It is assumed that, as long as it is proportional to, it is also proportional to. Therefore, an ideal can 
be estimated by approaching the number of cells to infinity. 
From the mesh convergence analysis, 660,528 cells were used to combine accuracy and minimize 
the computational costs. Figure 9 shows the aircraft pressure coefficient contours at the cruise conditions 
using this mesh. A surrogate model using the RMTB method was created for drag based on the SU2 


analysis in order to evaluate the aircraft performance in SUAVE. 





Pressure_Coefficient 
-].2e+00 0.5 0 0.5 1 1.2e+00 


Figure 9. Pressure coefficient contours of lower (left) and upper (right) aircraft surfaces. (« = 0°, M = 0.85). 


From this figure, one can observe that strong shock waves along the aircraft span on both the upper 
and power surfaces of the BWB aircraft, which results in large compressibility drag. Such behavior 
has not been captured by the low-fidelity analysis. In the low-fidelity analysis, a semi-empirical 
formulation was used to calculate drag due to compressibility. Although semi-empirical methods 
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provide good accuracy without computational costs, they may have limitations when the analysis is 
done for an unconventional configuration. Consequently, the higher aircraft drag predicted by SU2 
increased the fuel required to complete the mission by 24.1%. Therefore, the difference in fuel efficiency 
between the SE*A-LR, and B777 became equal to 10.6%, i.e., the new BWB aircraft has higher fuel 
consumption that Boeing 777. Table 7 shows a modified comparison between the SE7A-LR aircraft and 
B777 weights and fuel efficiency while using the higher-fidelity analysis. Appendix A presents details 
of the mission analysis using the results of the higher fidelity aerodynamic analysis. 


Table 7. Comparison of key aircraft parameters of reference long range aircraft SE*A-LR and B777 
using high-fidelity analysis. 





Parameter SE? A-LR (SUAVE-SU2) B777 Relative Change (%) 
MTOW (kg) 255,859.0 347,452.0 —26.4 
OWE (kg) 87,170.0 145,150.0 =599 
L/Dmax 43.0 21.0 104.8 
Cruise average L/D 5.0 18.5 —71.0 
Block fuel (kg) 135,681.0 109,290.0 24.1 
Fuel Efficiency (kg/seat/100 km) 3.01 2ilZ 10.6 


5. Aerodynamic Shape Optimization 


To mitigate losses that are created by compressibility drag, a discrete adjoint aerodynamic shape 
optimization using SU2 was performed. For the present optimization problem, the SLSQP algorithm 
available through the Python interface Scipy [40] was used. Only a single-point optimization for 
an average cruise condition was considered. The objective function was to minimize the aircraft 
drag, subjected to geometric constraints for 10 selected sections along with the aircraft. A Free-From 
Deformation Box (FFD) with 15 x 11 x 2 (total of 330) control points was used to modify the geometry 
at every optimization iteration. Table 8 shows the optimization problem statement and Figure 10 
shows the geometric representation of the FFD and constrained sections. 


Table 8. Objective function definition. 





Function/Variable Description Quantity 
minimize Cp Drag coefficient 
With respect to Z FFD control point z-coordinate 330 
Subject to Cr = 0.15 Lift coefficient constraint 1 
E> O:8tiacp Max thickness constraint for a given section 10 
Rye > 0.8RE,,<. Leading-edge radius constraint for a given section 10 





Figure 10. Free-Form Deformation (FFD) Box (red) and control sections (blue) of the SE*A-LR 
aircraft geometry. 


The lift coefficient constraint is based on an average lift coefficient during the cruise. The thicknesses 
of the initially selected airfoils satisfied the required volume for both passenger and fuel, however, 
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it was found that such a choice was conservative. The reduction of maximum thickness by 20% still 
satisfies volume requirements for fuel and passengers (although it negatively influenced the structural 
weight). Therefore, the thickness constraints for the optimization were defined to avoid getting a 
thickness reduction of more than 20% when compared to the initial design. Finally, a constraint of 
leading-edge radii reduction not exceeding 20% was imposed to avoid the possibilities of sharp-edge 
sections generation. 

Figure 11 shows a comparison between the initial and optimized designs. From this figure, one can 
observe that the shock wave strength was substantially reduced, which caused a reduction in inviscid 
drag by 77%. However, because the optimization has thickness constraints, shock waves have not been 
completely removed, so the compressibility drag is not equal to zero. Figure 12 shows a comparison of 
lift distributions between initial and final designs and elliptical lift distribution. The results show that 
the optimized geometry approached elliptical lift distribution when compared to the initial design 
which has a significant portion of lift produced by the fuselage and a mid-span loaded wing. However, 
from the structural perspective, elliptical lift distribution does not give the benefits from the structural 
point of view. Qin et al. [41] show that an elliptical/triangular lift distribution provides the best 
combination of minimum drag achievements and structural weight benefits. Present optimization has 
not captured this trend due to the absence of structural optimization coupling. Coupled aero-structural 
optimization will be performed during the next design iteration. 


Initial Design: Optimized Design: 
C, = 0.15 C, = 0.15 
Cy = 186.3 drag counts (Inviseid) | Cy = 42.5 drag counts (Inviscid) 
eg = —1.6752" ~ = —0.619° 











132 3 


Pressure Coefficient 
-].2e+00 0.5 0 0.5 1 1.2e+00 


10 


Figure 11. Pressure coefficient contours of the upper aircraft surfaces for initial and optimized configurations. 


Comparison of lift distribution 
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Percent span 


Figure 12. Comparison of lift distribution between initial and optimized designs. 


Figure 13 demonstrates the distribution of pressure coefficient along 10 control sections of the 
geometry. The blue lines show initial sections, while the red lines illustrate the modifications after 
optimization. From this figure, the smoothening of all the pressure distributions can be observed. 
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Besides, three trends can be observed: first of all, the maximum thickness of the inboard sections has 
not been substantially reduced, although a 20% reduction of thickness was allowed for all sections. 
The airfoil maximum thickness moved towards the leading edge making the aft part thinner. On the 
other hand, more significant thickness reduction is observed towards the wingtip. Second, the 
optimized wing demonstrates a 1 deg washout and increase the flight angle-of attack from —1.68 deg 
to —0.82 deg to reduce the shock wave strength, but still satisfying the lift coefficient requirements. 
Finally, a minor increase in the wing dihedral is observed. 
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Figure 13. Pressure coefficient and geometry sections. Blue curves represent the initial design and red 
curves represent the optimized design. 


6. Validation of the Drag Estimation Method using the Mixed Fidelity Model 


Because the Euler method is based on an inviscid flow assumption, its results lack the viscous 
drag. Although the lower fidelity method based on the flat plate analogy was used to compute the 
friction drag and add it to the inviscid drag of the Euler analysis, the accuracy of such a multi-fidelity 
method needs to be checked. For this purpose, a series of RANS analysis of the optimized geometry at 
the cruise Mach number of 0.85 at 35,000 ft for a fully turbulent wing was performed at different angles 
of attack. A mixed mesh that combines quadrilateral cells at the boundary layer and pyramid cells 
everywhere else was created while using Pointwise [42] software package. Figure 14 shows the mesh 
of half of the airplane. The mesh has seven-million cells, while the boundary layer was sized for Y+ 
of 1. SU2 was used for the Navier-Stokes equations solutions with the Spallart-Almaras turbulence 
model. Appendix A gives a summary of the pressure coefficient comparison between the Euler and 
the RANS solutions. Figure 15 shows a comparison of aerodynamic data between the RANS and the 
combined Euler and semi-empirical solutions. 
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Figure 14. CFD mesh for the RANS simulation. 
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Figure 15. Aerodynamic performance comparison between RANS and Euler with the flat plate analogy 
for the skin friction drag. 


The results show multiple trends. First, the RANS lift curve has a lower lift magnitude at a similar 
angle of attack when compared to the Euler solution. That is the cause of differences in pressure 
distribution and estimation of the shock wave strength. A similar trend was observed in [43] for the lift 
curve, where the lift coefficients at low angles of attack were similar with 0.5 deg difference between 
the experimental data and Euler solutions. After the angle of attack of 1 deg, the RANS lift curve 
starts showing non-linear behavior due to the generation of strong shock waves on the wing suction 
side. Such a trend is also observed for Euler solutions, but the Euler solution also does not simulate 
shock-induced separation, which happens in the RANS case; therefore, the reduction in the lift is lower 
in Euler compared to RANS.The drag polar shows that at low lift coefficient, the drag coefficient of 
RANS matches well with the Euler solutions plus the friction drag from flat plate analogy. However, 
deviations after the lift coefficient of 0.2 starts appearing. There is a good match between Euler (and 
friction drag) and RANS solution in drag polar until the lift coefficient of 0.5, but the drag polar of 
RANS then rises more rapidly due to the shock-induced separation. A similar increase in drag happens 
for lift coefficients below -0.1. The range of the whole mission lift coefficients falls between 0.1 and 
0.35, while the cruise segment is completed at lift coefficients between 0.12 and 0.2. For such a range, a 
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combination of Euler results with the flat plate analogy shows sufficient accuracy for capturing the 
aircraft aerodynamics and estimating its mission performance. 


7. Mission Performance Analysis of an Optimized Geometry 


After a single-point optimization, the mission was calculated again to obtain modified results 
with the account of high-fidelity aerodynamics. Appendix A shows this simulation results, while 
Table 9 provides a comparison of relative changes among low/medium fidelity, high-fidelity analysis of 
the initial aircraft and optimized high-fidelity results with respect to B777 results 


Table 9. Comparison of relative changes in key aircraft parameters with respect to B777. 


Initial SE*A-LR Low Initial SE7A-LR Optimized SE*A-LR 


ronemeren Fidelity (%) High Fidelity (%) High Fidelity (%) 
MTOW (kg) —58.5 —26.4 —51.1 
OWE (kg) 44,2 ~39.9 ~39.2 
L/Dmax 127.0 104.8 109.5 
Average cruise L/D 83.8 —71.0 21.6 
Block fuel (kg) —73.6 24.1 —55.5 
Fuel Efficiency (kg/seat/100 km) —76.5 10.6 —60.3 


The results show substantial improvements in mission analysis. However, improvements that are 
caused by new technologies do not match the low-fidelity analysis. Fuel efficiency improvement does 
not exceed 61% for the analysis while using higher fidelity tools. In addition, maximum L/D reduced 
from 50.5 to 39.5, which also demonstrates an increase in aircraft total drag. This shows that higher 
fidelity analysis plays an important role in initial stages of design and it should not be avoided if the 
design considers unconventional configurations. 

Although the results demonstrate substantial improvement in aircraft performances, there are still 
uncertainties in the results, which will be addressed in future works. Limitations of analysis are driven 
by the uncertainty of expected technological benefits and their real performance. It was assumed that 
80% of the wing will be laminar. However, no additional power required for the suction system was 
considered in the current analysis. Moreover, the effects of surface irregularities, such as landing gear 
doors, access panels, windows, etc., have not been considered. They might substantially reduce the 
laminarization of the vehicle or require unacceptable power for the suction system. In addition, the 
suction system weight was not considered at the moment due to a lack of information regarding the 
power required to suck the boundary layer from the wing. Other technologies, such as active load 
alleviation, boundary layer ingestion, and advanced structure concepts, also lack any surrogate data, 
which will be developed later under this research program. 


8. Conclusions 


The present manuscript describes a multi-fidelity design approach to design a fuel-efficient BWB 
aircraft. Multiple novel technologies have been integrated at the conceptual design phase to increase 
the aircraft efficiency and minimize negative environmental impact: active flow control (boundary 
layer suction), active load alleviation, boundary layer ingestion, and new materials and structure 
concepts. The presented approach was used to design a long-range BWB aircraft. 

A multi-fidelity approach has been used to increase design accuracy. The SUAVE aircraft design 
environment with the integration of novel technologies module was used for mission analysis, OpenVSP 
and Gmesh were used for geometry modeling and mesh generation, and AVL and SU2 were used for 
aerodynamic analysis. Euler equations CFD analysis was used to minimize the computational costs 
and capture transonic aerodynamic effects. 

A low/medium fidelity approach (using AVL for CFD) has shown significant improvement in 
fuel efficiency of a potential long-range aircraft—almost 76% increase in fuel efficiency as compared 
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to B777. However, mission analysis using SU2 demonstrated a substantial increase in drag due to 
compressibility that has not been captured by the lower fidelity analysis. A large increase in drag 
reduced fuel efficiency by 10.6% as compared to B777. 

High-fidelity discrete-adjoint aerodynamic shape optimization using SU2 was performed to 
minimize the adverse effects of the wave drag. FFD approach was used to modify the aircraft geometry, 
and an SLSOP algorithm was used for optimization. The objective function was to reduce drag for 
an average cruise condition and geometric constraints of a maximum thickness not being less than 
80% when compared to the initial design. The results showed a reduction of inviscid drag by 77% as 
compared to the initial configuration. 

The validation of the combination of Euler surrogate data and the flat plate analogy was also 
performed and the results used in SUAVE were compared to RANS simulations for cruise conditions. 
The results show sufficient accuracy between two methods for low lift coefficients. In addition, the 
cruise segment is performed at similar low angles of attack, so the method is sufficiently accurate to 
simulate the flight mission. 

A modified mission analysis with an optimized aircraft and higher fidelity aerodynamic analysis 
significantly improved the fuel efficiency as compared to the initial design: the fuel efficiency changed 
from +10.6% to —59.2% when compared to B777. However, an optimized solution did not match 
the lower fidelity analysis, which indicates the desire to introduce higher fidelity tools as early as 
possible in the conceptual design stage to have more accuracy during the design of unconventional 
aircraft configurations. 

Future work will focus on multiple tasks. First, a new design iteration with focus on stability and 
control, weights and balance, and system layouts will be performed. Second, the refinement of drag 
coefficient estimation needs to be performed: RANS analysis with embedded porous wall boundary 
condition responsible for the boundary layer suction will be performed, and new surrogate drag values 
will be obtained and used in SUAVE to more accurately determine potential fuel savings. In addition, 
surrogate models of structural weigh reduction, boundary layer ingestion, and active load alleviation 
will be implemented in SUAVE to more accurately estimate the potential benefits of novel technologies. 
Finally, an improved technique for the FFD definition will be implemented in SU2 to more uniformly 
distribute control points along with the aircraft geometry and increase the accuracy and flexibility of 
the shape optimization. 

The present manuscript only covers very preliminary results based on the methods developed 
in the first stage of the project. The major concerns include the actual capabilities of investigated 
technologies. Changes in technology assumptions may significantly affect aircraft performance and 
fuel savings. The design will face multiple iterations with progressively increasing levels of design 
and analysis accuracy, fidelity, and depth. 


Author Contributions: Conceptualization, 5.K., Y.L., A.E.; methodology, 5.K., Y.L., A.E.; software, S.K., Y.L.; 
validation, S.K.; formal analysis, 5.K.; investigation, S.K.; writing—original draft preparation, S.K., Y.L.; 
writing—review and editing, A.E.; visualization, S.K.; supervision, A.E.; funding acquisition, A.E. All authors 
have read and agreed to the published version of the manuscript. 


Funding: We would like to acknowledge the funding by the Deutsche Forschungsgemeinschaft (DFG, German 
Research Foundation) under Germany’s Excellence Strategy—EXC 2163/1-Sustainable and Energy Efficient 
Aviation—Project-ID 390881007. Furthermore, we acknowledge support by the Open Access Publication Funds of 
the Technische Universitat Braunschweig. 


Conflicts of Interest: The authors declare no conflict of interest. 


Aerospace 2020, 7, 87 


Appendix A Summary of SUAVE Analysis for All Analyzed Cases 


Table Al. Weight break-down of B777 calculated using SUAVE. 


Mass Items Unit Value 


Maximum Take-off Weight ke 347,452.0 


(MTOW) 
Operating Weight Empty (OWE) ke 157,432 
Wing kg 53,959 
Fuselage kg 29,032 
Empennage kg 7896.0 
Propulsion group keg 28,530 
Landing gear ke 13,898 
Systems (including opt. and furn.) 

: Breakdown of aie mass Kg oar 

Control systems kg 2691 

APU ke ers, 
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Figure A1. Aerodynamic performance of B777. 





17 of 27 


14000 


Aerospace 2020, 7, 87 


10 


Altitude (km) 


0 2000 


340000 


320000 


300000 


280000 


260000 


Weight (kg) 


240000 


220000 


0 2000 


Thrust (KN) 


Throttle 


= ° ° ° 
N os ron) © 


18 of 27 


S 


oo 
°o 
°o 


700 


fo] 
fo] 
o 


Ground Speed (km/h) 
ra 
oO 





0 2000 4000 6000 8000 10000 12000 14000 


550 
4000 6000 8000 10000 12000 14000 500 


450 


400 


Equivalent Airspeed (km/h) 


0 2000 4000 6000 8000 10000 12000 14000 





4000 6000 8000 10000 12000 14000 0 2000 4000 6000 8000 10000 12000 14000 
Range (km) Range (km) 
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Abstract: Morphing structures suitable for unmanned aerial vehicles (UAVs) have been investigated 
for several years. This paper presents a novel lightweight, morphing concept based on the exploitation 
of the “lever effect” of a bistable composite plate that can be integrated in an UAV horizontal tail. 
Flight dynamics equations are solved in Simulink environment, thus being able to simulate and 
compare different flight conditions with conventional and bistable command surfaces. Subsequently, 
bistable plates are built by using composite materials, paying particular attention to dimensions, 
asymmetric stacking sequence and total thickness needed to achieve bistability. NACAO011 airfoil is 
chosen for proving this concept. Wind tunnel tests demonstrate that the discrete surface is capable of 
withstanding the aerodynamic pressure. A remotely piloted vehicle is employed to test the discrete 
horizontal tail command during the take-off. The results show that, choosing a proper configuration 
of constraints, stacking sequence and aspect ratio for the bistable laminate, it is possible to tailor 
the snap-through mechanism. The proposed concept appears lighter and increases aerodynamic 
efficiency when compared to conventional UAV command surfaces. 


Keywords: morphing surfaces; bistable materials; RC-UAV 


1. Introduction 


Unmanned aerial vehicles (UAVs) may operate with various degrees of autonomy: either under 
remote control (RC) by a human operator or autonomously by on-board computers [1-3]. Compared 
to manned aircraft, UAVs were originally used for missions dull, dirty or dangerous and they can 
operate in environments contaminated by chemical, biological or radioactive agents. While they 
originally fitted mostly military applications, now their use is rapidly expanding to commercial, 
scientific, agricultural, recreational and other civil applications [4—8]. For this reason, in recent years 
there has been an increasing interest for (i) lighter-than-air UAVs with significantly higher endurance 
that can be used for persistent area surveillance [9]; (ii) long-endurance UAVs that can fly for several 
days [10-12] and (iii) small UAVs that are versatile, portable and easy to maintain [13,14]. 

In addition to these applications, UAVs provide an ideal platform for exploring and developing 
innovative solutions in the aeronautical engineering since the risk level is acceptable at this scale, 
and they have more manageable fabrication and less operational costs [15]. One example of this is 
the emergent idea of a morphing concept applied to aircraft lifting surfaces, which is the focus of 
this work. On this topic, several alternatives are reported in the technical literature, in terms not 
only of materials and structural parameters but also of feasibility and efficiency of the proposed 
solutions [16]. De Breuker and Werter in [17] demonstrated the strong influence of the aircraft speed 
on the maximum values of actuation forces in a morphing mechanism as well as the importance of 
the order in which various morphing deformations are executed. An adaptive trailing edge concept 
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was presented in [18] by combining conventional parts and innovative concepts, as the morphing skin 
assembled into a new structural system: an outboard wing trailing edge, consisting of three articulated 
(finger-like) rib covered by a multi—material hyperelastic skin. In [19] the antagonistic shape memory 
alloys (SMA) were used for wing morphing of small RC UAVs reducing the weight due to the actuation 
systems. Pecora et al. in [20,21] studied the design assessment of an innovative flap architecture for 
a variable-camber trailing edge: the reference geometry was based on a wing, where the conventional 
flap component was replaced by a morphing trailing edge based on compliant ribs. 

Among the morphing solutions proposed in last years, in terms of materials able to withstand 
aerodynamic loads, multistable materials were investigated for UAV morphing concepts: multistability 
is the property of varying geometrical shape with only small energy input and structures having this 
attribute represent an interesting candidate for the development of morphing solutions. Nicassio et al. 
in [22] explored bistable plates dynamic behavior for designing novel morphing command surfaces. 
In order to improve the UAVs aerodynamic performances, in [23] a prototype of morphing winglet 
utilizing a metallic bistable structure and polymer was presented. The metallic bistable structure was 
manufactured by bending a metal sheet and had two stable shapes: one was flat and the other was 
curved. To obtain an airfoil shape, the metallic bistable structure was placed in a mold with the shape 
of the airfoil and the polymer was poured into the mold for keeping the shape. The main finding of [24] 
is that multistable winglets are effective (especially at large dynamic pressure and angles of attack 
(AOAs)) in controlling the lateral and directional stability, but their efficiency in the lateral control, 
especially at low dynamic pressures, strongly diminishes respect to conventional ailerons. In addition, 
since their operational capabilities are strongly affected by the AOA, they can suffer roll reversal 
phenomena when wing’s AOAs are negative. Wang in [25] demonstrated that a morphing bistable 
composite structure could be produced, based on the principles of viscoelastically generated prestress: 
a composite plate structure was created in order to snap into one of two cylindrical states. In [26] 
cuidelines were proposed to design a morphing unsymmetric panel/actuator assembly, while in [15] 
a passive mobile surface with a bistable plate was presented. In the latter paper, the bistable morphing 
surface was activated using a proper configuration of constraints. These ad-hoc boundary conditions 
allowed the snap-though with activating forces of magnitude comparable with the differential pressure 
around a typical aeronautical airfoil. The possible applications of such behavior were the bistable 
composite integration in a low energy passive mobile surface able to autonomously react to pressure 
variations with a lift reduction when a maximum altitude was reached and vice versa. 

In the present manuscript a novel horizontal stabilizer concept for small RC-UAV (illustrated 
in the schematic of Figure 1) is designed by using composite bistable plates on which, in the future, 
micro fiber composite (MFC) will be applied and used as actuators. The main goals of this work are 
to (i) enhance the aerodynamic efficiency of the UAV control surfaces through the activation system 
integrated in the wing/tail box, (ii) replace the conventional UAV mechanical flight control system with 
an electronic one, with a resulting weight reduction, (iii) reduce the energy/load requirement since 
this system does not require transducers or servo-actuators to maintain the stable shapes but just to 


activate the snap-through. 
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Figure 1. Schematic of the present research. 
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The main advantages of this proposed concept are the following: the concept based on the 
activation of bistable is mechanically much simpler than the conventional mechanisms needed for 
the control of the mobile surfaces, consisting in command lines, pulleys, hinges and so on; one direct 
consequence of the simplicity is that also the inspection and maintenance are relatively easier for the 
bistable command system than the conventional one; in addition, the energy/load requested by the 
proposed solution is also lower than the one requested by the common kinematic chains, since in the 
latter the deflection of the mobile surface must be kept by an external reaction force able to withstand 
the aerodynamic pressure while in the former the deflection of the mobile surface is kept by an internal 
reaction force due to the thermal residual stresses inside the bistables. 

The present manuscript is not intended to cover all the design aspects typical of an aeronautical 
part even simple like the control surface of a small UAV. The activities hereinafter reported are 
a preliminary investigation about the technical feasibility of a control surface based on the employment 
of bistable materials. As such, the following research questions are addressed and answered: is it 
possible to control the flight of a typical UAV relying on a control surface capable only of discrete 
displacements? Can a bistable plate mimicking the conventional control surface with a reasonable 
weight and cost be designed and manufactured? Can such a plate actually be integrated into an 
aeronautical base profile? 

Following the choice of the UAV benchmark, the flight dynamics equations are implemented 
in MATLAB (Version R2015, The MathWorks, Inc., Natick, MA, USA) and Simulink environment in 
order to simulate different flight conditions (e.g., open and closed loop simulation, gust response and 
take-off maneuver). After checking the stability of closed loop bistable control system, the analytical 
results, in terms of discrete displacement of mobile surfaces to reach specific mission goals (i.e., elevator, 
ailerons and rudder), are used for the finite element (FE) simulations as design requirements (how 
much the bistable plate must deform). By using appropriate boundary conditions, it is demonstrated 
that the bistable shapes can be tailored to match these requirements. From the results of the FE analyses, 
carbon fiber reinforced polymers (CFRPs) plates were manufactured with aspect ratio (AR), total 
thickness, dimensions and asymmetric stacking sequence previously simulated in order to achieve the 
bistability. In the next step, wind tunnel tests are carried out in order to prove that the bistable system 
(integrated in a wing sample with the same UAV airfoil) remains in one of two stable geometries 
during the test and both the states are able to withstand the aerodynamic loads without any additional 
holding forces or locking mechanisms. Eventually, the discrete mobile surface behavior is validated by 
real flight tests. At this stage of the investigation using a conventional mobile surface in the “discrete 
mode” (just two possible deflections) is the only possible experiment. Following the experiments, 
the numerical model is furthermore validated matching the experimental take-off lengths with the 
numerical ones. The main finding is that the proposed innovative concept of bistable control surfaces 
actually works and the research can go to the next step: the implementation of an actual prototype of 
new horizontal stabilizer, by using bistable plates for the vehicle subject of study. 


2. Models, Materials and Methods 


2.1. Flight Dynamics Equations and UAV Features 


Under the hypotheses of flight in calm air, the motion equations valid for the UAV subject of 
investigation in the present study [27] can be expressed in an inertial frame and projected over a set of 
body-axes (see Figure 2). 
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Figure 2. Unmanned aerial vehicle (UAV) body axes (xp,yg,Zp) with: (u,v,w) velocities of the UAV centre 


of mass; (X,Y,Z) aerodynamic forces; (p,q,r) angular speeds; (L, M, N) roll, pitch and yaw aerodynamic 
moments, (7,9,W) roll, pitch and yaw angles; T thrust; Rg ground reaction; m mass. 


They form a system of 12 first order nonlinear ordinary differential equations (ODEs): 


e Dynamics of the centre of mass: 


u+qw—rv = (X+Tx)/m-gsind 
v+ru-pw = (Y +Ty)/m + gcos @sing ; (1) 
w+ po-—qu = (Z+Tz)/m-+ gcos@cosd—Rg 


e Attitude dynamics: 
ep — Inz(t + pq) - (ly - Lz )gr = £ 
yg -Ix2(r? -p?) - (le -Ie)rp = M_ (2) 
[zt — Ixz(p — qr) — (Ix —Iy)pq = N 
e Attitude kinematics: 
o =p+gqsingtan@-+rcos¢ tan 0 
@ = qcos@—rsing (3) 
my = qsingd/ cos@ + rcos@/ cos 8 


e ‘Trajectory of the centre of mass: 


Ry = ucos @ cosy) + (sind sin 6 cos) — cos @ sin) + w(cos sin @ cos W + sin sin ) 
Re = ucos Osiny + v(sing sin @ sin) + cos cos w) + w(cos@ sin @ sin — sing cos w) 
h = —(-usin@ + vsingcos 6 + weos@cos @). 
(4) 


In order to carry out the experimental campaign, a homemade RC UAV was selected. The choice 
of this UAV was based on the following reasons: it is simple to be managed; its basic main structure 
allows the future use of bistable surfaces; its aerodynamic behavior is easy to be simulated and many 
investigations on this type of aircraft are reported in the literature and research works; lastly the 
installed airfoil is widely used in mini UAV and small UAV, that represent the final application of this 
study. Geometric dimensions of the RC UAV are reported in Figure 3, mass, aerodynamic and thrust 
features in Table 1. 

In order to verify the centre of mass position and to obtain the inertia moments of the actual UAV, 
a CAD model was developed using CATIA software (Version V5-6 R24 R2014, Dassault Systemes SE, 
Velizy-Villacoublay, France). 
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Figure 3. UAV geometric drawing, dimensions in millimeters. 


Table 1. Unmanned aerial vehicle (UAV) main features. 








Specification Symbol Value Specification Value 
Mass and Inertia Aerodynamics 
Mass M 6.5 kg Wing and Tail Profile NACA 0011 
L 0.681 kg m2 Thrust 
ly 1.375 kg m? DLE Engine Power 5 CV 
2 . . 

iehy snaincaio lacie I, 2.0121 kg a DLE Engine Displacement 30 cc 
ley —0.001 kg m Propeller Diameter 0.48 m 
les 0.022 ke m? Propeller Pitch 0.30 m 
ly 0.0001 kg m2 - - 


To evaluate the aerodynamic characteristics of the NACAO0011 airfoil (wings, horizontal and 
vertical tail), XPFOIL software [28] was used. This is an interactive program widely used for the design 
and analysis of subsonic isolated airfoils. The main airfoil characteristics are reported in Figure 4 and, 
according to the low speed considered in this work, Reynolds and Mach numbers were found to be 
equal to 400,000 and 0.04 respectively for all simulations. 

The data in Figure 3, Table 1 and Figure 4 were used to perform in Simulink the following studies: 


e Open-loop simulations by applying a specific command input to the elevator, ailerons and rudder: 


the command consisted in a 1 degree positive step; 


e Closed-loop simulations (gust response by using Dryden wind turbulence model) with a full-state 
feedback control by means of the linear quadratic regulator (LOR) method; 


e ‘Take-off length simulations. 


The last two simulations were carried out with both continuous and discrete operation modes of 
the mobile surface, to compare their effect on the UAV motion. 
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Figure 4. XFOIL aerodynamic properties of NACAO011. 


2.2. FE Simulations and Manufacturing of Bistable Plates 


A thermomechanical FE model [29] was developed to design the bistable composite plate as 
a morphing surface able to provide the same deflections of the continuous mobile surfaces. The entire 
analyses were split up in two load cases: the thermal simulation to induce the curved shapes and the 
transient simulation of the snap-through. The commercial software ANSYS® (Transient Structural 
Module, version 2019 R3, Ansys Inc., Canonsburg, PA, USA) was used [30]. To create the bistable 
stabilizer (for replacing the conventional one in Figure 2), a bistable rectangular plate with an AR 
of 2 (80 mm X 160 mm) was simulated by stacking four laminae following a (02/902) sequence and 
exploiting the “lever effect” [15]: the leverage amplifies the input force to provide a bigger output effect. 
The ratio of the output force to the input is the mechanical advantage of the lever. This mechanism 
allows us to gain mechanical advantages by moving the septum, in terms of plate displacements 
and activation forces: to increase the displacement and decrease the activation load, the distance 
between the constrained edge and the septum must be as small as possible (see Figure 5a). The actual 
characteristic boundary conditions were carried out with one constrained edge and the adjustable 
septum (as in Figure 5b). 








/ FE plate 
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Figure 5. (a) Lever effect graphic and (b) finite element (FE) model (two stable states) on a UAV stabilizer. 
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The FE model consisted of shell elements. Convergence studies on the mesh size were carried out 
in order to achieve an acceptable accuracy of the results. The option “Large Deflection” was employed 
in each step because large deformations are associated with small displacement increments of the plate 
during the process (geometric non-linearity). 

The available material was the prepregs Hexply 8552 with unidirectional carbon fibers IM7 (used 
in primary aerospace structures [31]). The main mechanical properties shown in Table 2 filled the 
engineering data module of the software for all simulations. 


Table 2. Hexply 8552 data. 


Technical Data Value 
Nominal cured ply thickness 0.131 mm 
Nominal laminate density 1570 kg m~? 
Young’s modulus 0° direction 165 GPa 
Young’s modulus 90° direction 11.4 GPa 
Coefficient of thermal expansion 0° direction 6x10-7C"! 


Coefficient of thermal expansion 90° direction 2.86 x 10° C7! 


The numerical results provided by FE models allowed us to know: (i) the displacements evaluated 
at the free end of the plate for different septum positions, (ii) the maximum deformations and (iii) the 
activation forces. The displacements values allowed us to evaluate the deflection capabilities of the 
discrete mobile surface. The knowledge of plate deformations was useful since it drives the choice of 
possible locations for a proper use of MFC P1/F1 actuators [32]: these actuators, in fact, bonded onto 
the plate and electrically activated, are deformed inducing the bistable snap-through. The activation 
forces intensity is strictly related to the amount of charge requested by the MFC. 

In order to obtain experimentally the discrete shapes of bistable plates, four laminates with four 
layers of Hexply prepregs (see Figure 6) were manufactured. 





(a) | (b) 


Figure 6. Bistable laminates manufacturing: (a) Hexply 8552 prepregs on the mold and (b) vacuum bag. 


The samples were cured following the conventional cycle suggested by the manufacturer, setting 
the temperature at 180 °C, with the overall curing cycle of 3h. In order to increase the thermal residual 
stresses and to obtain a higher level of curvature in the samples, the laminates, after the curing was 
completed in the autoclave, were immediately brought out so that a fast cooling down occurred in 
open air at room temperature. After the cooling down, each laminate was observed to have the two 
expected stable states (in Figure 7). 





Figure 7. Cont. 
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(b) 


Figure 7. Experimental bistable plates: (a) first stable state and (b) second stable state. 


3. Results and Discussion 


3.1. Simulink Results 


In this subsection the UAV response to input commands considered in the previous section 
was analyzed. The open-loop simulations were run to demonstrate that the Simulink model was 
effective. If a unit step is assigned to a mobile surface, its effect on the flight dynamics will be the 
same regardless if the command will be continuous or discrete. The closed-loop simulations instead 
provided information on how the control actions from the controller were dependent on the system 
dynamics. The closed-loop simulations were performed to make a comparison between continuous 


and discrete bistable command surfaces providing information useful to understand the technical 
feasibility of a discrete command. 


3.1.1. Open-Loop Simulations 


Unit step inputs after 1s of trimmed flight were considered in all cases as inputs to the control 
mobile surfaces (initial trim conditions corresponding to an horizontal flight with speed and altitude 
equal to 10 m/s and 50 m respectively), in order to validate the Simulink model by means of flight 
mechanics equations. 

Figure 8 represents the dynamic response in terms of motion variables to a positive (pitch-down) 
step input on the elevator. The first 20s of simulation are reported, where it is possible to recognize 
the short-term response of short period variables, namely a and g. The AOA rapidly reaches a new 
equilibrium with well damped oscillations, whereas the pitch rate converges towards a small value. 
When a positive elevator step is applied, the initial force increment coming into play is directed 
upwards, and it causes a small but strictly positive increment to the flight-path angle. This causes the 
aircraft to enter an accelerating dive and u increasing. 
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Figure 8. UAV response to the unit step elevator. 
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The response to a positive aileron step command was characterized by an initial transient 
(5-6 s) dominated by the roll mode, during which the UAV achieved a negative roll rate (Figure 9). 
The increment in the roll angle achieved during this phase caused the aircraft to sideslip (6 became 
negative). A steady state yaw rate r was achieved for a constant roll angle ~, which means that the 
aircraft entered a steady turn. From the practical standpoint, it should be noted that the aileron control 
power in this case was related to the ailerons placed close to the root of the half-wing in order to limit 
the excitation of wing bending and torsional modes. This explains the small value of the roll rate p 


achieved in the maneuver. 














































































































“0 5 10 15 20 25 30 0 10 20 
time (s) time (s) 


Figure 9. UAV response to unit step ailerons. 


When a positive step command was applied to the rudder, the first effect was a linear acceleration 
around the yaw axis and the growth of the sideslip angle 6, which became positive (Figure 10). 
The sideslip angle generated a negative rolling moment, thus inducing an increment of the roll angle 
~, as achieved by the aileron command, but related to a different source. The response was clearly 
dominated by sizeable Dutch-—roll oscillations, superimposed over a slow spiral convergence. At the 
steady state, also in this case, the UAV was significantly banked. 
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Figure 10. UAV response to the unit step rudder. 
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3.1.2. Closed-Loop Simulations 


Discrete displacements of mobile surfaces give the possibilities of performing a wide range of 
maneuvers. Simulink was used to develop a stability and control augmentation system suitable for the 
chosen UAV: the comparison between continuous and discrete commands on closed-loop performance 
was made. As in [33], a LOR was developed for the virtual model, with full-state feedback on: 


e (apes — a) and (qpgs — g), with apgs and gpgs equal to arryy and to zero respectively, for 
longitudinal stability; 
e § (ppes — p) and (rpges — r), with ppgs and rpgs equal to zero for lateral—directional stability. 


The subscript DES refers to the target value. The resulting closed-loop behavior for a 3D gust 
disturbance (duration 30 s in Figure 11) on the desired value of «, g, p and r is presented by the 
following pictures. 











Uc (m/s) 








We (m/s) 























Wo (m/s) 


“0 5 10 15 20 25 30 
time (s) 


Figure 11. Dryden wind turbulence velocities. 


By using a Simulink function, it was possible to implement the Dryden wind turbulence model 
(continuous) block that uses the Dryden spectral representation to make turbulence acting on the UAV 
model (in Figure 11 the component velocities of the used Dryden gust were plotted). As it is also 
possible to see from the figure, the speed components were not correlated with each other and never 
exceeded +5 m/s. Since the turbulence in a real scenario was highly dependent on the considered 
altitude, for the specific case subject of study, the pertinent turbulence model was adopted considering 
that the evolution of the UAV was enclosed in an altitude range between 45 and 55 m. Another aspect 
was the turbulence duration: most of the actual gusts lasted a few seconds; in the case of the present 
investigation a gust lasting about 30 s was used for the simulations representing the worst dynamic 
scenario. The gust components were incorporated directly into the nonlinear dynamics equations 
simulating UAV behavior, like perturbations of the aircraft velocity. 

The LQ regulator was thus synthesized for the model: the cost functions, previously defined, 
were associated with the deviations due to the gust perturbation of UAV altitude and attitude from the 
trim values. The control law developed by LOR obviously does not take into account if the command 
is continuous or discrete, but simply provides a deflection value of the surface needed for achieving 
the DES (i.e., desired) conditions. The resulting closed-loop gust response with continuous mobile 
surfaces is provided by the dotted lines in Figure 12a. In order to satisfy the cost functions, the resulting 
elevator deflection falls in the range between +20° (Figure 12a top); the ailerons instead show very 
low deflections needed to guarantee the desired flight conditions (Figure 12a centre); the rudder must 
deflect +10° to maintain the lateral directional stability (Figure 12a bottom). The closed-loop behavior 
of the conventional LOR controller was acceptable, with minor changes around the trim conditions 
(these results are not reported in Figure 12b for the sake of clearness of the plot). These results confirm 
the good stability performance of the chosen UAV. 

A control law was then developed for the UAV model operated by discrete commands: +2° for 
the elevator and +0.5° for the ailerons were the discrete limit deflections. The resulting control law 
is represented by continuous blue lines in Figure 12a where the deflections were evaluated in order 
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to satisfy the cost functions. In this case the closed-loop behavior was also clearly acceptable (see 
Figure 12b): the UAV trajectory almost always remained within a 5 m diameter cylinder, centered at 
the trim trajectory. More in details, in Figure 12b the different trajectory components were reported: 
z represents the altitude and (x,y) plane was the horizontal plane parallel to the ground. From the 
Figure 12b it is clear that the vehicle, under the action of the gust and the reaction of the discrete 
controls, moved around the trim conditions with apparent oscillations that were always confined in 
about 5 m. To simulate also the actual bistable devices, the piezoelectric MFC actuator was integrated 
into the Simulink model taking into account its activation law. For the snap transition between the two 
stable positions of the control plate, it is necessary to send an electric discharge to the MFC actuator. 
Therefore, a charging time, allegedly affecting the overall dynamics of the system, was added to the 
control law. The charging time was needed by the condenser for the subsequent discharge generation 
activating the MFC: the numerical simulations including this charging time (order of 1s) demonstrated 
that this aspect was not critical for the control of the vehicle. 

The results showed in this section proved that the discrete controller concept was actually working. 
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Figure 12. UAV gust response: (a) elevator, ailerons and rudder deflections (continuous and discrete 
behavior) and (b) UAV 3D trajectory with discrete command surfaces. 


3.1.3. FEM vs. Experimental Plate Results 


The numerical FE model was developed (i) to numerically evaluate the plate geometry after 
the cooling down, (ii) to simulate the plate displacements in order to compare these results with the 
experiments and (iii) to obtain the maximum cambers with different degrees of application of the “lever 
effect”. FE models furthermore drive the choice of the bistable dimensions and stacking sequence 
to replace UAV conventional mobile surfaces. A plate 80 mm x 160 mm was simulated by stacking 
four laminas (02/902) of Hexply 8552 (Hexcel, Stamford, CT, USA). Typical results (deformations) 
are shown in Figure 13: the numerical maximum cambers are also reported and compared with the 
experimental results. 

The “lever effect” was introduced to constrain the laminate used like an aerodynamic mobile 
surface (see Figure 14): “S” indicates the distance in mm from the septum position to the bistable 
constrained edge (in order to better understand the meaning of the symbol “S” used throughout the 
text, Figure 5 can be looked at). At each value of “S” it is possible to evaluate the maximum deflection 
of the discrete mobile surface (vertical displacement f, rotation angle). 

The airfoil NACAO011, currently used for the stabilizer of the chosen UAV, was virtually integrated 
with the design bistable plate: the discrete deflections of the bistable plate (a) and of the compliant 
profile (b) are plotted in Figure 14. Table 3 summarizes the geometric results: these proved that by 
choosing proper configuration of constraints, lay-up and AR for the bistable plates it is possible to 
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tailor the snap-through mechanism to make bistables actually operating like mobile surfaces. For 
each S value, the bistable discrete mobile surface had a deflection angle (in Table 3) comparable with 
a conventional UAV control surface. 
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Figure 13. FE bistable model (two stable state shapes deformations): simulated and measured cambers. 
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Figure 14. Discrete mobile surface with “lever effect”: (a) plate deflection f with different value of S; 
and (b) airfoil NACA0011 with different bistable angles 6 and deflections f. 


Table 3. Deflection f and 6 for each values of S. 


S (mm) Deflection f (mm) Deflection Angle 6 (°) 





40 84.1 ~28 
80 56.9 ~20 
120 2/ ~10 


3.2. Wind Tunnel Experiment 


The FE simulations showed that, using the “lever effect”, the deflection of an aerodynamic surface 
made of bistable can be comparable with the one of conventional configuration. In order to prove 
experimentally this idea, a wing sample was made up of balsa, using NACAO011 as a base profile. 
A bistable plate exploiting the lever effect was integrated in the base profile (Figure 15a). The forward 
spar was used as clamp for the short side of the plate, while the intermediate support was realized by 
a mobile septum in order to achieve the lever effect before presented and discussed. In this way the 
MFC actuator in future developments of this concept will be able to activate the snap-through: each 
stable state of the composite laminate provided a discrete position of the mobile control surface, due to 
the compliant skin at the profile trailing edge (Figure 15b). 
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(a) (b) 
Figure 15. Experimental NACA 0011 airfoil (a) sample without the skin and (b) discrete mobile surface 
OFF (top) and ON (bottom). 


The capability of the wing sample to keep the shapes at several flight speeds and at different 
AOAs was evaluated by carrying out specific tests represented in Figure 16. 





(b) 
Figure 16. Wind tunnel experiment: (a) 1. Sample profile with discrete mobile surface, 2. angles of 


attack (AOA) tool (min —15°, max +20°), 3. pitot tube, 4. test chamber and 5. wind speed regulator and 
(b) negative (LEFT), neutral (CENTRE) and positive (RIGHT) deflection of discrete mobile surface. 


The scope of the aerodynamic campaign was to demonstrate the stability of every state of the 
mobile surface when operating in an external field. In this context stability of an initial state means 
that, following a small perturbation of the external flow (i.e., caused by a local separation of the flow 
around the profile or by the turbulence), the bistable keeps its geometric shape without snapping to the 
other shape. The possible bistable initial states were four, two for each stable state (see Figure 17, where 
it is clear that the plate can be installed in one specific geometric configuration and the corresponding 
reverse one, UP and DOWN respectively). During the experimental tests the external field was 
accelerated at different AOAs waiting for a “natural” snap through, i.e., not activated by any sensor 
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(MFC). This campaign was extremely useful for understanding if the mobile surface can be prone to 
unwanted deflection. 


Stable state 2 






Stable state 2 








Figure 17. Bistable plates configurations. 


An accurate aerodynamic experimental campaign was carried out in a wind tunnel of the type 
“open circuit”. S was equal to 80 mm, to which corresponds a deflection angle of about 20°, as it can be 
seen in Table 3. Two different sample configurations were tested, according to the scheme in Figure 17, 
to allow the study (varying AOA and flight speed) of the snap behavior, due to the aerodynamic load, 
for positive and negative deflections of the mobile surface. Table 4 lists the speeds at which the snap 
through occurred: the stable state 1 (in both the configurations) was the most suitable to withstand the 
aerodynamic external pressure, whereas the stable state 2—UP configuration suffered negative AOAs 
and the stable state 2—DOWN configuration suffered positive AOAs, ie., snapping without external 
activation. This result was consistent with the static characterization of the bistable plates reported 
in [15] where it is experimentally shown that the stable state 1 is stiffer (i.e., requires more energy for 
the transition to the other state) than stable state 2. 


Table 4. Experimental wind tunnel results. 


Initial Stable State 
AOA (°°?) 1—UP 2—UP 1—DOWN 2—DOWN 








Snap Speed (km/h) 

+20 85 - - 35 
LO - - - 45 
+10 - - - 50 
+5 - - - 65 

0 - 60 - 70 
—5 - 70 - - 
—10 - 55 - - 
—15 - 45 - - 


In conclusion, the proposed morphing concept passed the wind tunnel test, showing significant 
aerodynamic withstanding at AOAs and flight speeds that are typical of the size of the employed UAV. 
These results also suggested a solution to replace the conventional UAV stabilizer (with weight around 
0.5 kg) with an innovative horizontal tail with discrete bistable elevators (see CAD in Figure 18a). 
The proposed concept appeared lighter (bistable plates and lever effect tool weigh around 0.3 kg) 
and increased aerodynamic efficiency (without any servomechanism outside of the tail airfoil) when 
compared to conventional UAV command surfaces. 
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Four bistable plates (in Figure 18a) could be installed on a modified leading edge of the new 
horizontal tail with the lever effect tool integrated (Figure 18b): a pair of plates could be installed 
in the UP configuration, and a pair of plates could be installed in the DOWN configuration (see 
Figure 17). In this way when they are not activated, their first stable shape provides the neutral 
position of the stabilizer; when one pair is activated it is possible to achieve a deflection of the stabilizer 
upwards/downwards depending on which one of the two pairs (UP/DOWN respectively) is activated. 





(a) (b) 


Figure 18. New horizontal tail concept: (a) bistable horizontal tail in CATIA and (b) UAV equipped 
with the concept of a horizontal bistable stabilizer. 


3.3. Takeoff Length: Numerical vs. Experimental Results 


This final subsection summarizes the results of three experimental tests and Simulink simulations 
related to the UAV take-off length. The experimental campaign was carried out with the conventional 
UAV stabilizer operating in “discrete” mode (just two possible deflections): this choice of using the 
convectional stabilizer and not the bistable one was motivated by the circumstance that, in the event 
of failure or an unwanted dynamic evolution of the UAV, the conventional commands would have 
allowed a prompt and full recovery of the UAV. As suggested by the RC UAV manufacturer, during 
the maneuver and the Simulink simulations, the engine throttle underwent a gradual increase up to 
60% of the maximum in 10 s for the take-off maneuver. In the following pictures, white-red reference 
elements were placed at a distance of 5 m from each other. 

The first take-off test was carried out assuming the stabilizer deflection equal to —10°; this 
deflection was constant from the beginning of the test until the UAV took off from the ground: the 
vehicle was soared in flight after having covered 18.7 m of runway in 4 s (Figure 19). 





Figure 19. Runway length at first takeoff time. 


The second take-off test was carried out assuming the stabilizer deflection equal to 0°; this 
deflection was constant: the UAV was not able to take-off using the whole runway length of 35 m. 
The last take-off test was performed like a traditional maneuver: taxing with the stabilizer in neutral 
position until 3 s, and subsequent —10° deflection on the horizontal tail command. Under these 
conditions, the UAV was able to lift off after 12.5 m of runway, in 3.5 s. 

Table 5 reports the numerical results obtained Simulink in which the motion equations were 
solved and the experimental results. 
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Table 5. Take-off results. 








Test Results Experimental Simulink Error (%) 
ar aS Runway (m) 18.7 17.6 5.9 
7 a Time (s) 4 3.8 5.0 
=e Runway (m) 35 35 - 
& ery Time (s) 6 5.7 5.0 
3 0-3 s, dp = 0° Runway (m) 12.5 13.6 —8.8 
from 3s on, dg = —10° Time (s) 3.5 Sad Oud 


The errors on time and space were so explained: 


e Ground friction coefficient: in Simulink simulations this was constant, in actual tests the runway 
was definitely uneven and not properly described by a constant friction coefficient; 

e Measurement time: the error of human reaction was about three tenths of a second; 

e Wind: the actual weather was quickly changing, with small gusts from different directions, while 
no gust model was included in the Simulink code; 

e UAV weight distribution: it was assumed homogeneous in the Simulink model for sake of 
simplicity, and at this stage of the research this small error can be accepted; 

e Ground effect: it was not included in the Simulink model. 


In conclusion, the Simulink code was able to predict the UAV behavior, either during taxing and 
flight maneuver, either with conventional and discrete mobile control surfaces. Now it can be stated 
that the proposed innovative concept of bistable morphing control surfaces on a RC-UAV can go to the 
next research step: implementation of the actual prototype of new UAV horizontal stabilizer, by using 
composite bistable plates integrated with the lever effect tool due to the data gathered during this 
research phase. 


4. Conclusions 


In this study a new horizontal stabilizer was presented for RC-UAV application, by using composite 
bistable laminates. The nonlinear flight dynamics equations were implemented in Simulink to study 
different flight conditions. These numerical results were used as inputs to the FE model of bistable 
composite plates. By using the “lever effect”, it was demonstrated that the shapes of the bistable 
plates could be tailored to mimic the conventional command surfaces behavior. Composite bistable 
plates were built with AR, dimensions, asymmetric stacking sequence and total thickness, previously 
simulated. The new bistable elevator was tested in the wind tunnel and the discrete mobile surface 
behavior was validated through real flight tests, carried out using a actual RC-UAV. The present 
manuscript just reports a preliminary investigation on the technical feasibility of a control surface based 
on the employment of bistable materials. Basic research questions aimed at proving that the concept 
could actually work were addressed positively. Future work will be based on a robust design of the 
control surface consisting in the optimization of the aerodynamic efficiency of the surface, through 
numerical and wind-tunnel tests, and in a proper evaluation of the energy required for the operation of 
the control. After all the design stages will be completed, a prototype of bistable horizontal stabilizer, 
activated by a MFC actuator, will be installed on a UAV and flight tested. 
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Abstract: In order to solve the accidents caused by aircraft icing, electro-impulse de-icing technology 
was studied through numerical simulation and experimental verification. In addition, this paper 
analyzed in detail the influence of the number, placement arrangement, and starting time of pulse 
coils on the de-icing effect, which plays a guidance role in the design and installation of the 
subsequent electro-impulse de-icing system. In an artificial climate chamber, the new de-icing 
criteria were obtained by tensile test, and the platform for the electro-impulse de-icing system was 
built. Replacing the skin with an aluminum plate, an electro-impulse de-icing system with a single 
coil was used. A three-dimensional skin-ice layer model was established by using Solidworks 
software. The finite element method was adapted. Through comparison between the de-icing 
prediction results and the test results in the natural environment, it was proven that the calculation 
process of de-icing prediction was correct, which laid a theoretical foundation for the selection of 
the number, placement arrangement, and starting time of the pulse coils. Finally, in this paper, 
by choosing the leading edge of NACA0012 wing as the research object, the influence of the number, 
placement arrangement, and starting time of pulse coils on the de-icing effect was analyzed. The results 
show that to get a better de-icing effect, the electro-impulse de-icing system with two impulse coils 
should be selected. The two coils were installed in the central position of the top and bottom surfaces 
of the leading edge, respectively. In addition, one of the impulse coils started working 1200 us later 
than the other one. 


Keywords: aircraft; wing; system; ice; icing; deicing; ice mechanics; electro-impulse; coils; 
structural dynamics 


1. Introduction 


Since the earliest days of aeronautics, icing was found to be a crucial problem for aircraft flight. 
In-flight ice accretion occurs on the leading edge of an aircraft wing and usually covers only 2% of the 
wing chord, with the thickness of the ice layer being about a few centimeters. However, even an ice 
layer of a few centimeters thickness at the key parts of the aircraft is enough to cause flow separation 
and destroy lift, increase drag and reduce the maximum lifting capability, affect the control surface 
effectiveness, and in some cases decrease engine performance and stability [1-3]. Therefore, people have 
been working on aircraft anti-icing methods for many years [4-6]. The present ice protection methods 
are hot bleed air, freezing point depressants, and electro-thermal resistance heating. However, they all 
have potential limitations, for example, pneumatic boots bonded to the ice prone surface are subject to 
corrosion and damage by external objects, therefore, they need to be replaced every two or three years. 
The freezing-point depressant systems achieve the purpose of aircraft anti-icing by releasing ethylene 
glycol through many orifices in the leading edge of an aircraft. There are two main hazards: on the 
one hand, it increases the weight of the aircraft; on the other hand, since ethylene glycol is a toxic 
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substance, release of this substance into the air can cause environmental pollution. The principle of 
electro-thermal systems is to use resistance pads to heat the aircraft to above the melting temperature 
of ice, and in general it requires dedicated generators, which results in significant cost. 

The aforementioned methods cannot meet the performance requirements for new aircraft. 
Therefore, it is imperative to develop a safe and reliable de-icing system. A new method for 
electro-impulse de-icing has emerged. 

Electro-impulse de-icing (EIDI) is one of the mechanical de-icing systems, and it assures the safety 
of aircrafts in an icing condition. It has major advantages, such as low energy, minimal maintenance, 
ereat reliability, and low cost and weight [7,8]. 

The basic circuit of the electro-impulse de-icing system is shown in Figure 1. The working 
principle is the pulse coils are connected to a high voltage capacitor by low resistance, low inductance 
cables. When the switch is turned on, the discharge of the capacitor through the impulse coils creates a 
rapidly forming and collapsing electro-magnetic field. According to Maxwell’s law, we know that the 
time-dependent magnetic field induces eddy currents in the metal skin. Therefore, the instantaneous 
impulse force of several hundred pounds magnitude is obtained by the Lorentz force formula, but the 
duration is only a few hundred microseconds. A small amplitude, high acceleration movement of the 
skin acts to shatter, de-bond, and expel the ice. 








SWITCH 


ENERGY STORAGE 
CAPACTTOR 
SKIN 


Figure 1. Electro-impulse system basic circuit. 


Researchers [9-11] have studied the phenomenon of ice layer peeling on the surface of the skin 
under the action of a single coil system. The number of coils is greater than one when a real EIDI 
system is installed on the aircraft [12-14], so the position arrangement and start-up time of the multiple 
coils may affect the vibration condition to the effect of the ice removing effect. The de-icing criteria 
used in previous studies [15-17] are from the literature [18]. The adhesion between the ice layer and 
the skin is different under different icing conditions. Therefore, the conclusions of [18] are not suitable 
for all icing environments. Ice tensile strength and ice/aluminum adhesive shear strength used in [19] 
are 1.5 Mpa and 0.5 Mpa, respectively. However, this does not explain how the values are obtained. 
Therefore, the study of de-icing criteria is very important for the electro-impulse de-icing system. 

All the experiments in this paper were obtained in an artificial climate chamber. Replacing the skin 
with an aluminum plate, the electro-impulse de-icing system with a single coil was used. The de-icing 
range obtained by experiments after each pulse was compared with the calculated result, and it was 
proven that the calculation process for de-icing prediction was correct. In addition, in the same icing 
environment, the adhesion between the ice layer and the skin was measured by a tensile test, and the 
new de-icing criteria was obtained. Finally, a three-dimensional NACA0012 wing-ice layer model 
was established; the effect of the number, placement arrangement, and the start-up time of pulse coils 
on the de-icing ratio was analyzed, which provided a reference for the design and installation of the 
subsequent electro-impulse de-icing system. 
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2. System Composition and Structural Dynamic Studies 


The structural dynamics equation is expressed as [20] 
IMJU + [C]U+ [K]U =F, (1) 


where [M] is the overall quality matrix, [C] is the overall damping matrix, [K] is the overall stiffness 
matrix, F is the impulse force, and U is the node displacement, U is the node accelerated velocity, U is 
the node velocity. 

In the continuous linear elastic body vibration, the relationship between normal stress, shear stress 
and displacement and external force is described. The expression for the stress can be written as 


[o] = [D][B}{uy’. (2) 


Substituting the node displacement into Equation (2), 0x, Oy, Oz, Txy, Tyz, Txz for each node can 
be obtained. 

The characteristic equation for the stress state is solved by the theory of one-dimensional 
cubic equation. Therefore, the practical calculation formulae for the principal stress can be written as 





i, F (cos 0 =) 
07 = 3 3 S35 ool } 
I - 0 
03 = 3 F (cos os + V3sin3 5) 
3Ip_I? QI, In -2 PP -2713 Oe 
where P = —,+; q = ——z7—; 9 = arccos -4(-5) (0<O0<7). 


From the formulae above, 01, 02,03 can be obtained. Substituting 01, 02, 03 into Equation (3), 
the equivalent stress can be written as 


g = [4 1(o1 - 09)? + (09-03)? + (03-01) 3) 


3. Calculation Method Verification 


3.1. Geometric Model 


In this paper, an aluminum plate was used as the research object, and the material was the 
same as the wing skin. The research of electro-impulse de-icing systems mainly includes three parts: 
electrodynamics research, structural dynamics research, and de-icing prediction. The electro-dynamics 
studies of this EIDI included the calculation of impulse current, the study of the magnetic field behavior, 
and the calculation of the impulse force of each grid point on the surface of the aluminum plate, 
which provided the basis for structural dynamic research. Structural dynamics studies include stress 
calculations, which provide the basis for de-icing prediction. The model was established by using 
Solidworks, with a length of 420 mm, width of 420 mm, and thickness of 1.5mm. The pulse coil was 
located on the lower surface of the aluminum plate, and the gap between the impulse coil and the 
skin was 1 mm. It had an inner diameter of 20 mm and an outer diameter of 128 mm. The material 
characteristics are shown in Table 1. 

The structural dynamics associated with electro-impulse de-icing have proven to be a difficult 
and challenging problem. The structural dynamics study is mainly to calculate the equivalent stress 
between the ice layer and aluminum plate under the action of pulse force. If the total pulse force is 
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applied to the aluminum plate, the calculation accuracy will be affected. Therefore, the aluminum plate 
was divided into different regions in the radial direction. Then, the impulse pressure at different radial 
points and different times was placed on the corresponding regions on the surface of the aluminum 
plate. This greatly improved the accuracy of the de-icing prediction. It is assumed that the ice layer 
with 4 mm thickness was evenly covered on the aluminum plate with 2.0 mm thickness. The material 
properties are shown in Table 2. 


Table 1. Material characteristics. 








Material Relative Permeability Conductivity (S/m) Density (Kg/m°) 
Aluminum plate 1.000021 1.74 x 10” 2780 
Impulse coil 0.999991 5.8 x 107 8933 


Table 2. Material properties. 





Material E/GPa V o/(kg-m?) 
Aluminum 7.1 0.33 2700 
Ice 5.5 0.3 897 


The size of the meshing affects the calculation results of structural dynamics. In this paper, 
the modal analysis of the aluminum plate was calculated by the finite element method. The calculated 
result of first-order natural frequency was 46.5 Hz. Considering the boundary conditions and bending 
vibration of the aluminum plate, the aluminum plate was regarded as a free Bernoulli-Euer beam 
model at both ends, and its natural frequency theoretical calculation formula [17,18] is as follows: 


Tt Le [EI 
fn = ("+ 5) ry ai (4) 


where | is the length of the aluminum plate, A is the cross-sectional area, p is the density, I is the 
moment of inertia of the section, and E is the elasticity modulus. 

The first-order natural frequency obtained by the finite element method was 46.5 Hz, and the 
theoretical calculation was 44.5 Hz. The relative error between the two was very small, only 4.3%. 
It indicates that the meshing sizes above were feasible, which provided a basis for the correct structural 
dynamics calculation. 


3.2. De-Icing Criteria 


The key factor of de-icing prediction research is the de-icing criteria. According to [9], there are 
two main de-icing criteria. One is to consider only a single strength factor at the ice-skin interface, 
and the other is to consider the tensile and shear stresses at the ice-skin interface. However, no matter 
which kind of the de-icing criteria is used, the equivalent stress and the shear stress are both determined 
values (oy = 1.44 MPa, ty = 0.4 MPa). This is contrary to the fact that the equivalent stress and the 
shear stress between the ice layer and aluminum plate is different under different icing environments. 

In this paper, the equivalent stress and the shear stress between the ice and aluminum plate was 
measured by a tensile test. Figures 2 and 3 show the icing box and tension test device, respectively. 
The icing test was carried out in an artificial climate chamber. As shown in Figure 4, the height and the 
diameter of the artificial climate chamber were 11.6 m and 7.8 m, respectively. The wind speed in the 
built-in wind tunnel can be up to 10 m/s. The sprinkler system consisted of 14 high-pressure water 
mist nozzles. 
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Figure 3. Tension sensor. 





Figure 4. Artificial icing test. 


The temperature curve during whole test was —20 degrees. Note that if a large amount of water is 
to be poured into the icing box at first, the water easily overflows. This causes icing around the outer 
surface of the icing box, resulting in inaccurate measurements. Therefore, a better option is to pour a 
little water into the icing box at first until the water in the icing box freezes completely, and repeat this 
procedure until the icing box is full of ice. In the natural icing environment, the outer surface of the 
icing box can also freeze. Before the tensile test, the icing around the icing box should be carefully 
removed so that the measurement result is accurate. 

In addition, to improve measurement accuracy, the aluminum plate should be polished with gauze 
to remove the oxide film and rust on its surface. During the measurement process, the operating speed 
was controlled so that it accurately read the sensor readings. The results for the adhesion between the 
ice and the aluminum plate are shown in Table 3. 

The tangential force is produced by electro-impulse de-icing system is small, therefore, this paper 
only needed to consider that the ice layer falls off when the tensile stress between the ice and aluminum 
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plate is greater than the maximum tensile stress. As shown in Table 3, the maximum tensile stress 
Omax between the ice and aluminum plate is 


EF 
ax = — = 0.15 MPa, 


where Frave is the average value of the adhesion and S is the area of the ice box, 20.25 cm?2. 


Table 3. The adhesion between the ice and the aluminum plate. 














Test Object Equivalent Stress (kN) 
(1) 0.323 
(2) 0.325 
(3) 0.274 
(4) 0.299 
(5) 0.289 
(6) 0.252 
(7) 0.310 


The maximum stress Omax between the ice and aluminum plate obtained in this paper and the 
result obtained by [21] are not much different. Therefore, it is indicative that the maximum stress Omax 
measured by tensile test is correct. 


3.3. Result Analysis 


3.3.1. Analysis of Experimental Results 


The icing test was carried out in natural icing station. The icing results of aluminum plates are 
shown in Figure 5. Icing thickness was 4 mm, as shown in Figure 6. 

The circuit parameters are: voltage: UO = 780 V; charging capacitor: C = 990 uf; inductance: 
L=19uH, Lw = 5 HH. The de-icing results of skin are shown in Figure 7. 

It can be seen from Figure 7 that the de-icing rate reached 73.6% after three-pulse. It is indicative 
that the electro-impulse de-icing system was feasible. 





Figure 5. Icing results of the aluminum plate. 
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Figure 7. Ice layer shedding process. 
3.3.2. Analysis of Calculation Results 


A three-dimensional impulse coil-aluminum plate model was established by using Solidworks 
software. The importing of the model into the Maxwell software is shown in Figure 8a. An external 


circuit shown in Figure 8b was applied at the cross section of the impulse coil. 


ModelName 





(a) 3D impulse coil-aluminum plate model. (b) An external circuit. 


Figure 8. Electro-dynamic simulation model. 


Under the action of the instantaneous impulse force F, the equivalent stress between the skin and 
the ice layer can be obtained by Equations (1)-(3). The de-icing range on the surface of skin was found 
by using the new de-icing criteria. In order to improve the calculation result accuracy of de-icing 
prediction, the structural dynamic response was obtained by using the finite element method. The force 
at different times in different positions of the skin was applied to the corresponding position, as shown 
in Figure 9. The tangential magnetic field B, and the eddy current density Jeqay at different positions 
and at different times of the aluminum plate were obtained by simulating. Therefore, the instantaneous 
impulse force can be obtained by Equation (5), as shown in Figure 10. 


F = { J eddy X Bdv, (5) 
V 
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where F is the impulse force, Jeqay is the eddy current density, and B is the magnetic field. 








Figure 9. Excitation and loading area. 


FZz/N 


300 400 
t/ps 





Figure 10. Relationship between force and time. 


As shown in Figure 10, the distribution of force increased first and then decreased along the radial 
distance. The time required to reach the peak pressure was 100 us, which meets the requirements of 
narrow pulse width. 

According to the process of calculation, the equivalent stress was obtained by the finite element 
method, as shown in Figures 11 and 12. 

The de-icing ratio obtained by the experiment was 73.6% after three pulses, and the de-icing ratio 
obtained by calculation reached 72.1% after three pulses. Although the de-icing ratio was very close, 
there was a difference in the de-icing range. The experimental results show that the lower and left side 
of the ice layer on the surface of the aluminum plate were completely detached, and the calculation 
results did not fall off. This was due to the effect of the external factors, such as wind speed. During the 
experiment, the ice on the top of the skin was easily able to fall off under the action of external wind 
speed. This factor was not considered in the calculation process of de-icing prediction. 

In summary, the de-icing ratio obtained by the experiment was 73.6% after three pulses, and the 
de-icing ratio obtained by calculation reached 72.1% of the de-icing prediction after three pulses. In the 
range of errors permitted, the comparison between the experimental results and the predicted results 
indicated that the calculation process of de-icing prediction was correct. 
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Figure 11. Equivalent stress. 





Removed ice 


De-icing 72.1% 


Figure 12. De-icing prediction results. 


3.3.3. The Effect of Ice Thickness on the De-Icing Ratio 


To study the changes of de-icing ratio with ice thickness while maintaining all other parameters 
unchanged, we assumed that the thicknesses of ice is 1 mm, 6 mm, 8 mm, and 10 mm, respectively. 
The de-icing ratios obtained by calculation are shown in Table 4. 


Table 4. Effect of ice thickness on de-icing ratio. 


Ice Thickness (mm) 1 6 8 10 
De-icing ratio 80.7% 71.04% 53.44% 10.42% 


It can be seen from Table 4 that the thicker the ice is, the smaller the de-icing ratio will be while 
maintaining all other parameters unchanged. 


4. Example 


The number of pulse coils was greater than one when a real EIDI system was installed on the 
aircraft, so the number, position arrangement, and start-up time of the multiple coils may have affected 
the vibration condition to the effect of the ice removing effect. The NACA0012 wing with a chord 
length of 1000 mm was taken as the research object in this paper. Since the wing was covered with 
ice mostly in the leading edge part [22-24], it was only necessary to analyze the de-icing ratio of the 
leading edge of the wing. A three-dimensional wing leading edge-ice layer model was established by 
using Solidworks, as shown in Figure 13. Its length along the chord was 300 mm, the length along the 
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exhibition was 600 mm, and the span of one bay of the leading edge was 300 mm. The thickness of the 
leading edge was 1.5 mm and the thickness of the ice layer was 1.0 mm. 





(a) Electro-impulse de-icing system with four coils. (b) Electro-impulse de-icing system with two coils. 
Figure 13. Position arrangement of multiple pulse coils. 


4.1. Position Arrangement of the Multiple Coils 


In order to study the influence of multiple coils on the de-icing ratio, the single coil simulation 
method was extended to a multiple coil simulation. Regarding the possibilities for weight reduction of 
EIDI systems, it was also important to determine the minimum number of coils, that is, whether it 
is necessary to place a coil in every bay of the leading edge. Because the established model was 
symmetrical, it only needed to select two or four pulse coils. The position arrangements of the four 
and two pulse coils are as shown in Figure 13a,b, respectively. The structure dynamics simulation 
of the four coil system was carried out. If only two of four coils were used, it was divided into three 
kinds of modes, which are respectively Mode 1, Mode 2 and Mode 3, as shown in Table 5. Mode 4 is 
the two coil system. 


Table 5. Four kinds of two coils. 


Mode (1) (2) (3) (4) 
Position arrangement 1,20r3,4 1,30r2,4 1,40r2,3 5,6 


The de-icing results of simultaneous vibration of the four coils are shown in Figure 14. For the 
electro-impulse de-icing system of the four coils, the de-icing results of the simultaneous vibration are 
shown in Figures 15-17, in which only two of four coils were used. The de-icing results of simultaneous 
vibration of the two coils are shown in Figure 18. 
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Figure 14. Electro-impulse de-icing system of four coils. 
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Figure 16. 1,3 or 2, 4 of four coils are used (Mode 2). 


1.116 Max 
0.99203 @ 
0.86808 
0.74414 
0.6202 

0.49625 
0.37231 
0.24837 
0.12442 0.00 250.00 500.00 (mm) 
0.00047897 Min 125.00 375.00 


Figure 17. 1, 4 or 2,3 of four coils are used (Mode 3). 
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Figure 18. Electro-impulse de-icing (EIDI) system of two coils (Mode 4). 


It can be seen from Figures 14—18 that the de-icing ratios varied with the number or the position 
arrangement of pulse coils. In addition, the ice layer fell off within a certain distance from the central 
position of the pulse coil. This was because the magnetic field was stronger in this range. The stronger 
the magnetic field, the bigger its magnetic force, thus, the ice layer was easy to fall off. The de-icing 
ratios of the several models above of the electro-impulse de-icing system are shown in Table 6. 

It can be seen from Table 6 that for the electro-impulse de-icing system with two pulse coils, 
the de-icing ratio of Modes 1, 2, and 3 were only 36.19%, 38.35%, and 41.43%, respectively. The de-icing 
ratio of Mode 4 was 47.22%, which was the maximum. The de-icing ratio of the electro-impulse 
de-icing system with four coils was 69.83%. According to the experimental results, the ice layer at the 
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boundary of the skin would also fall off under the influence of external factors. Therefore, from the 
perspective of not affecting aircraft flight and low energy consumption, the position arrangement of 
Mode 4 is optimal for the research objects mentioned above. 


Table 6. De-icing ratio comparison. 


Electro-Impulse De-Icing Electro-Impulse De-Icing System with Two Coils 
System with Four Coils Mode 1 Mode 2 Mode 3 Mode 4 
De-icing ratio 69.83% 36.19% 38.35% 41.43% 47.22% 


4.2. The Effect of Start-Up Time of Pulse Coil 


In the case of Mode 4, the starting times of the pulse coils 5 and 6 were different, and the de-icing 
ratios obtained by the calculation were different. Take as an example the system with two pulses, 
pulse coil 5 started to work at 0 us, and pulse coil 6 started to vibrate from 0 us, 100 us, 600 us, 
and 1200 us, respectively. The vibration de-icing result at 0 us was the same as that of Mode 4. 
The results of the ice layer detachment of pulse coil 6 starting to work from 100 us is shown in 
Figure 19b. Figure 19c shows the pressure waveform at the coil radial direction r = 30 mm. The results 
of the ice layer detachment of pulse coil 6 starting to work from 600 us and 1200 us are shown 
in Figures 20 and 21, respectively. 
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Figure 19. De-icing results of pulse coil 6 starting 100 ps later than coil 5. 
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Figure 20. De-icing results of pulse coil 6 starting 600 us later than coil 5. 
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Figure 21. De-icing results of pulse coil 6 starting 1200 us later than coil 5. 


It can be seen from Figure 19 that the de-icing ratios obtained by calculation were different under 


the different starting times of pulse coil 6. De-icing prediction results are shown in Table 7. 


Table 7. De-icing ratios. 


Start Time Delay of Coil 6 with 
Respect of Start Time of Coil 5 


De-icing ratio 47 22% 51.38% 45.42% 54.03% 


0 us 100 us 600 us 1200 us 


It can be seen from Table 7 that if pulse coil 6 started working at the same time as coil 5, or 600 Ls 


later, the de-icing ratios obtained by calculation were almost the same. This was because the interval 


between the continuous applied force at the same position of the wing was the same. However, if pulse 
coil 6 started working 1200 us later than coil 5, the de-icing ratio was the largest, which served as a 


euiding function for the subsequent pulse coil vibration control. 


5. Conclusions 


(1) 


In an artificial climate chamber, the maximum adhesion between ice layer and aluminum plate was 
0.15 Mpa obtained by the tensile test. Since the tangential force generated by the electro-impulse 
de-icing system was very small, as long as the equivalent stress between the ice layer and the 
aluminum plate was greater than 0.15 Mpa, the ice layer fell off, which laid a foundation for the 
simulation calculation of the subsequent electro-impulse de-icing system. 

We replaced the skin with an aluminum plate and built a platform of the electro-impulse de-icing 
system in an artificial climate chamber. The range of de-icing of the skin surface obtained by 
the experiments was recorded under different pulse times. A three-dimensional aluminum 
plate-ice layer model was established by using Solidworks. The results obtained by the three-step 
calculation of electrodynamics, structural dynamics, and de-icing prediction were compared 
with the experiment results of the electro-impulse de-icing system. It was indicative that the 
calculation process of de-icing prediction obtained by this paper was correct. 

Taking the leading edge of the NACA0012 wing with a chord length of 1 m as the research object, 
its length along the chord was 300 mm, the length along the exhibition was 600 mm, and the span 
of one bay of the leading edge was 300 mm. The leading edge had a thickness of 2.0 mm and the 
ice thickness was 4.0 mm. The de-icing ratio of the four pulse coil system was 69.83%. Under the 
action of the two pulse coil system, the de-icing ratio of Mode 1, 2, and 3 was lower, and the 
de-icing ratio of Mode 4 reached 47.22%. From the perspective of low energy consumption and 
without affecting the safe flight of the aircraft, the aforementioned structure was best to select the 
electro-impulse de-icing system with two pulse coils. 

With regard to the electro-impulse de-icing system with two-pulse coils, the starting time of the 
pulse coil was different, and the de-icing effect was different. If pulse coil 6 started working at 
the same time as coil 5 or 600 us later, the de-icing ratio was almost the same. This was because 
the interval between the continuous applied pressure at the same position of the wing was the 
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same. However, if pulse coil 6 started working 1200 us later than coil 5, the de-icing ratio was the 
largest, which served as a guiding function for the subsequent pulse coil vibration control. 
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